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ABSTRACT 

A  set  of  formulas  proposed  by  Moritz  and  constituting  a  second  order  solution 
of  Molodensky's  problem  has  been  transformed  in  a  way  that  the  resulting  formulas 
are  numerically  tractable.  The  difficulty  in  the  original  formulas  arose  from  the 
occurence  of  singular  integrals.  Their  regularization  introduces  derivatives  of 
gravity  anomalies  and  terrain  heights  and  necessitates  a  numerical  differentiation 
procedure.  Spline  function  interpolation  has  been  chosen  to  deal  with  this.  The 
existing  truncation  theory  has  been  extended  to  cover  more  sophisticated  truncation 
procedures  as  well  as  a  succession  of  heavier  and  heavier  smoothed  versions  of  the 
integrand. 
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INTRODUCTION 


Moritz  (1969)  proposed  a  set  of  formulas  which,  under  certain  assumptions, 
constitute  a  second  order  solution  for  Molodensky's  problem. 

These  formulas  are  numerically  not  tractable  because  singular  integrals  and 
even  iterated  singular  integrals  occur.  The  occurence  of  the  singular  integrals  implies 
that  the  formulas  are  integro-differential  formulas  in  disguise.  The  singularities  can 
readily  be  removed  by  allowing  derivatives  to  show  up  openly.  Green's  second  formula 
for  surfaces  is  the  essential  tool  to  regularize  the  integrals,  and  regularization  should 
take  place  only  within  a  small  cap  centered  at  the  point  of  interest.  For  the  outer  zones 
the  original  form  of  the  integrals  is  mo  re  advantageous  since,  in  many  cases  the  kernel 
tapers  off  quickly  with  increasing  distance. 

For  remote  areas  a  smoothed  version  of  the  involved  functions  like  gravity 
anomalies  and  terrain  heights  can  be  used.  The  existing  truncation  theory  has  been 
reviewed  and  somewhat  extended  to  cover  more  sophisticated  truncation  procedures 
as  well  as  a  succession  of  heavier  and  heavier  smoothed  versions  of  the  involved 
functions.  Appendix  B,  which  is  concerned  with  truncation,  is  expected  to  be  of  some 
interest  beyond  the  immediate  applications  in  this  report. 

The  occurence  of  derivatives  in  the  modified  formulas  introduces  the  necessity 
of  numerical  differentiation.  Assuming  that  the  functions  are  only  given  at  discrete 
locations  or  in  block  average  form  the  spline  interpolation  approach  has  been  chosen. 

This  interpolation  method  has  some  advantages  over  the  classical  methods  of  poly¬ 
nomial  or  trigonometric  interpolation.  The  bi-cubic  spline  interpolation  procedure 
assumes  that  the  functions  are  given  at  discrete  locations  forming  a  rectangular  grid. 
These  can  be  accomplished  by  prediction  methods  applied  to  an  originally  irregular 
pattern  of  sample  points.  In  the  case  of  block  averages  the  most  points  of  the  blocks 
can  be  used  to  form  the  rectangular  grid 

This  work  has  been  carried  out  under  some  time  pressure.  Therefore  extensive 
numerical  tests  could  not  be  performed.  Some  tests  have  been  made  with  bi-cubic  spline 
interpolation  and  also  with  one  of  the  singular  integrals. 
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The  proposed  methods  are  certainly  time  consuming.  However  a  second  order 
approach  will  be  employed  if  in  a  limited  area  a  good  accuracy  is  needed.  In  this  case 
the  chosen  approach  appears  to  be  computationally  feasible. 
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THE  ORIGINAL  FORMULAS 


We  shall  deal  with  equ.  (8. 1)  to  (8.3)  in  Moritz  (1969).  In  somewhat  different 
notation  these  formulas  read  as  follows: 

a(§)  =  3WG  J  St<?’ "»?)  (Ag(r?)  +  G^rj)  +  G2(t? ) )  drXr))  - 

r 

2 

1  r  (h(n)-h(f))  a  _  / _ \  v  /i  1\ 


—  f 

4nRG  " 


(s-n) 


Ag(ti)dT(n) 


(1.1) 


v(^)  =  -  J*  Grad^St( § •  tj)  (Ag(Tj)  +  G1  (tj)  +G2(ti) )  dF^)  + 

r 

+  iSPc  / ' h(5) >2°rad?7^> ag(”,dr(’’)  ■ 


^(SJ+G^S) 


Grad  h(5) 


where 


G1<5)  *  4  J 

g2<5)  =  m  J  h<,.3'h<f  Gi<?i)dr<i)  + 
r  je  ( 5’r? ) 

+  ^TAg(5)  Grad  h(|)  2 


(1.4) 


Explanation  of  notation: 


unit  sphere 

unit  vectors,  denoting  points  on  T.  They  play  the  role 
of  spherical  coordinates: 

l  =  (sin  6  cos  X,  sin  0  sin  X ,  cos  9) 
astronomical  coordinates  are  obtained  by  putting 

_  n  a 

<P  ~T  9 

inner  product  of  %  ,  tj:  frequently  denoted  by  cos  ip 
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&  (5  •  v )  - 


St(S-Tf)  = 


R 

G 

Ag(?) 


h(5) 

a(5) 

v(5) 


Gradsf(  *)  = 


distance  between  5  and  rj: 

£(^■77)  =  £(cos  0)  =  2  sin 

Stokes’  function 

St(l-Tj)  =  St(cos  0 )  =  — L__  +  .  .  . 

sin  2 

(See  Appendix  B,  equ.(B.24)) 
mean  earth  radius 
mean  gravity  value 

gravity  anomaly  referring  to  a  point  on  the  earth’s 
surface  with  astronomical  coordinates  implied  by 

height  above  sea  level  in  f . 

height  anomaly  at  5  in  the  sense  of  Molodensky  et  ai 
(1962);  in  the  literature  usually  denoted  by  Also 
the  undulation  of  the  quasi -geoid. 

deflection  of  the  vertical;  viewed  as  a  surface  tan¬ 
gential  vector  to  F.  In  die  literature  usually  de¬ 
noted  by  its  components  * ,  77  in  a  localized  coord¬ 
inate  system. 

denotes  the  surface  gradient  of  a  function  f(^)  de¬ 
fined  on  L  The  surface  gradient  is  a  vector  tan¬ 
gential  to  T  in  5  and  pointing  into  the  direction 
of  maximal  increase  of  f(5).  The  length  of  the 
vector  equals  the  rate  of  the  increase. 


A  derivation  and  discussion  of  the  physical  meaning  of  these  formulas  will  not 
be  given  here.  They  will  be  used  as  a  starting  point  for  this  study  except  for  one 
thing:  During  the  derivation  of  the  formulas  "planar  approximation''  has  been  used. 
This  means  that  in  deriving  a  formula  for  a  certain  quantity  a  relative  error  of  the 
order  h/R  has  been  tolerated.  We  shall  do  the  same  whenever  there  will  be  a  need 
for  modifying  some  of  the  formulas . 

The  accuracy  which  shall  be  aimed  at  is  about  three  digits  in  the  quantities 
n(c)  anil  v(p).  This  does  not  necessarily  mean  that  this  accuracy  can  be  guaranteed 
ilirougliout.  It  means,  rather,  that  wo  do  not  worry  about  approximations  which 
create  errors  of  less  than  three  digits  in  a(£)  and  v(£). 
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2.  THE  INTEGRO- DIFFERENTIAL  CHARACTER  OF  THE  FORMULAS 

The  formulas  (1.1)  to  (1.4)  are  frequently  referred  to  as  integral  formulas. 
Though  it  is  recognized  that  some  of  the  involved  integrals  are  singular,  the  evaluation 
of  the  formulas  is  usually  viewed  as  a  problem  of  numerical  integration. 

Let  us  first  have  a  look  at  the  integral  (1. 3)  serving  the  evaluation  of  G  j(  §  ). 

This  integral  is  in  an  abbreviated  way  written  as 

Gi  -  mr  Jipr-Sidr  (2-‘> 

r 


Although  we  shall  use  for  the  actual  evaluation  of  this  integral  another  transformation; 
we  write  (2.1)  now  as 


Gi 


1  a  h  Ag  -  h  Ag 

2itR  J 

r 


dr 


(2.2) 


Here  we  have  two  integrals  of  the  type 


(2.3) 


The  spherical  harmonics  equivalent  of  this  is 

ynm  =  -  n  (2.4) 

This  shows  an  amplification  of  the  higher  harmonics  of  the  same  type  as  it  occurs  in 

a  differentiation.  As  outlined  in  Meissl  (1971)  a  transformation  of  the  type  (2.4)  can 

k  + 1  k 

be  viewed  as  a  transformation  from  a  space  H  into  a  space  H  where  the  menr 
k+ 1 

bers  of  H  possess  certain  generalized  derivatives  up  to  order  k  +  1  whereas 

Jr 

those  of  H  possess  only  such  derivatives  up  to  order  k. 

Thus  (2.3)  acts  similar  to  a  differentiation  procedure  which  always  makes  things 
more  rough  than  they  have  been  before. 


Daring  the  evaluation  of  G2  such  a  roughing  effect  occurs  twice.  Stokes' 

k  k  +  1 

formula  has  an  opposite  effect.  It  transforms  from  H  into  H  .  It,  so  to 
speak,  decreases  the  roughing  procedure  by  one  step.  Nevertheless  in  summary  it 
turns  out  that  a(§)  depends  in  some  way  on  derivatives  of  h  and  Ag.  Whereas  the 
mean  contribution  to  a(£);  namely, : 


R 

4nC 


J  St  Ag  d  I‘  , 

r 


is  smooth  (one  step  smoother  than  Ag),  the  correctional  terms  are  rougher  .  This 
has  serious  consequences  on  any  theoretical  discussion  of  the  Molodensky  approach 
which  will,  however,  not  be  undertaken  here. 

The  consequences  which  have  to  be  dealt  with  here  are  that  there  will  be  a 
need  to  evaluate  in  a  consistent  way  derivatives  of  functions  which  originally  are 
given  only  in  disci ete  locations.  We  shall  briefly  outline  now  what  we  mean  by  con¬ 
sistent.  Take  for  example  Green's  second  identity 


J (f  Lap  g  *  g  Lap  f )  d  ?  =  J*  f  (Grad  g,  v  )  -  g  (Grad  f,  v )  daB 
B  3  B  *** 

Lap  denotes  the  surface  Laplacean  operator  with  respect  to  F.  B  is  a  sub-area  of 
T  with  boundary  SB.  The  unit  vector  v  is  tangential  to  T  and  normal  to  a  B. 

(v  is  directed  outward  of  B). 

If  now  the  functions  f  and  g  are  given  only  at  certain  discrete  locations  then 
no  integral  in  (2.5)  can  be  evaluated.  If  the  involved  derivatives  are  computed  by 
some  crude  and  inaccurate  numerical  interpolation  and  differentiation  procedure 
then  (2.5)  could  possibly  not  be  verified.  We  shall,  however,  use  transformations  of 
our  integrals  which  are  based  on  (2.5).  Therefore  what  is  clearly  needed  is  a  pro¬ 
cedure  to  evaluate  the  derivatives  of  the  functions  in  a  way  that  (2.5)  holds  to  a  sat¬ 
isfactory  degree  of  accuracy.  We  have  to  interpolate  the  discrete  values  in  a  way 
that  the  function  is  continuously  differentiable  up  to  the  second  order.  There  are 
many  ways  to  do  this  and  the  result  is,  of  course,  not  unique.  We  reject  trigonometric 
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or  polynomial  interpolation  in  the  classical  sense  and  choose  instead  the  spline 
function  approach.  The  latter  yields  generally  a  smoother  surface  because  local 
disturbances  have  hardly  any  influences  on  distant  areas.  This  is  more  fully 
explained  in  Appendix  A.  The  way  data  are  assumed  to  be  given,  i.e.  what  the 
locations  of  the  discrete  arguments  are,  is  explained  in  section  4.  For  the  immed¬ 
iately  following  sections,  it  suffices  to  know  that  the  functions  and  their  derivatives 
up  to  second  order  are  defined  everywhere  on  r. 


3.  REGULARIZATION  OF  THE  INTEGRALS 

A  singular  integral  can  be  transformed  into  a  regular  one  if  the  integrand  is 
sufficiently  smooth.  This  may  for  the  first  be  illustrated  by  the  following  simple 
example . 

3.1.  A  one -dimensional  example 

Take: 

g(x)=j£&Zldy  (3.1) 

(x  -  y  r 

-l 

We  assume  that  <p(x ,  y )  is  twice  continuously  differentiable  and  that 

<p(x,  x)  =  0  (3.2) 


The  integral  (3.1)  does  not  exist  if  the  integrand  is  replaced  by  its  absolute  value. 
The  integral  is  therefore  singular.  Writing  instead  of  (3.1) 


x-e 

1 

+  1 

ge(x) 

*  1  +  J  - 

(e)  J  .  say. 

-1 

x-* 

-  1 

the  singular  integral 

g(x  )  is  defined  by 

lim  g£(x). 

€~*0 

Applying  partial  integration  toward  gf  (x  )  yields 


g  (x )  =  gfc  y)  €  +  fiM.) 

x  -  y  x  ~y 


Applying  partial  integration  a  second  time  gives 


+  1  /  v 

,  <Pv(x>y) 

(€)J  ^ - dy 


x 


- 1 


x-  y 


gcOO 


(p(x,  y) 

x-y 


x  -  e 

+ 

-  1 


<p  (x,y ) 

x  -y 


1 

x+e 


. . .  .(con't) 
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+  <py(x,  y)  £n|x  -  yl 


X -€ 


+  <py(x,  y )  £nfx  -  y| 


1 

x+c 


•  (O  f  <Pyy(x,y)  £nlx  -yl  dy 


-  1 


If  we  now  let  c—  o  then  the  various  poles  cancel  out  neatly  and  we  obtain 


g(x)  =  _  *  -  *  <x.-l)ln|l+x|  +  <py(x,l)  fn  |l-x| 

+  1 

'  J  <^yy(x»y)  £a>x  *  y 1  dy  (3-3> 

- 1 

The  integral  is  now  completely  regular,  which  means  that  it  exists  even  if  the  integrand 
is  replaced  by  its  absolute  value. 

3.2  Regularization  of  M-type  integrals 

Let's  go  now  to  the  unit  sphere  T.  Assume  an  integral  of  the  form 

g<5>  =  l  7^')dr<,,,  (3’4) 

with 

<p(5.5)  =  0  (3.5) 

and  </>( 5 ■  77 )  at  least  twice  continuously  differentiable.  We  call  this  a  type  M-inte- 
gral  since  the  first  order  correction  term  in  Molodensky's  formula  is  of  this  type,  e. g. 
equation  (1. 3).  The  singular  integral  is  defined  after  excluding  a  very  small  circular 
cap  of  half  opening  angle  e  around  §  and  then  letting  e  -*  o. 

The  equivalent  of  the  above  two-fold  partial  integration  is  in  the  two  (or  more) 
dimensional  case,  nothing  but  Green's  second  formula  (2.5). 

We  transform  the  integral  not  over  the  whole  sphere  F  but  only  over  a  certain 
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spherical  cap  C  surrounding  5  and  having  half-opening  angle 


It  is  shown  in  Appendix  C  that 


-  S  -  w  ■ 


f  {*(§• 


r?)  S.  (cos  ^0) 


- - -}  Lap,,  <p(l,v)dr(v)  " 

s  wJ  v 


f  <P(l>r))  daC(rj) 


cos  ~2~ 

2  .  - 
S.  (cos  0Q)  ac 

±  f  dr(v) 

4  £  4<§-»7) 


(3.6) 


The  formula  becomes  simpler  if  the  cap  is  extended  over  the  whole  sphere; 
in  that  case  the  line  integral  vanishes .  However  it  is  not  advisable  to  do  so .  The 

Q 

term  l/H  ( §  -  77)  in  (3.4)  tapers  off  quickly  for  larger  distances.  This  is  not  so  in 
(3 . 6)  where  only  terms  like  l/f(5*tj)  occur.  Therefore  for  distant  areas  (3.4) 
offers  greater  advantage.  The  term  Lap^fpd'Tj)  in  (3. 6)  certainly  tends  to  fluctuate. 
If  the  integral  is  extended  over  too  large  a  cap  then  the  positive  and  negative  contri¬ 
butions  of  lap^d-T})  certainly  cancel  to  some  extent.  Therefore  the  cap  should 
not  be  too  large  allowing  only  a  few  positive  and  negative  extremes  of  Lap^ 
within  its  domain. 


In  all  numerical  approaches  toward  the  singular  integrals  of  Physical  Geodesy 
which  are  known  to  me  the  integration  according  to  (3.4)  is  carried  out  very  close  to 
the  point  g  of  singularity.  For  a  very  small  remaining  cap  some  rather  crude  pro¬ 
cedure  Is  used  which  involves  some  kind  of  numerical  differentiation.  This  does  not 
seem  to  be  a  very  efficient  procedure.  A  singular  integral  yields  infinity  if  the  inte¬ 
grand  is  replaced  by  its  absolute  value.  Therefore  if  only  a  very  small  cap  around 
the  singularity  is  excluded,  there  are  necessarily  large  positive  and  negative  contri¬ 
butions  which  nearly  cancel.  This  is  numerically  disadvantageous  even  if  a  precise 


10 


analytical  expression  for  the  integrand  is  available.  An  analytical  representation  for 
the  integrand  appears  to  be  a  necessity  in  any  case.  If  one  is  available  then  derivatives 
may  be  formed.  Formula  (3.6)  can  then  be  evaluated  withe,,  t  any  numerical  trouble 
resulting  from  mutual  contributions  of  large  quantities  of  opposite  sign.  The  cap 
shall  not  be  too  small  but  also  not  teo  large  as  we  have  pointed  out  earlier. 

M-type  integrals  appear  in  our  formulas  of  section  1  multiplied  by  a  factor 
i.  In  that  case  the  last  integral  in  (3.6)  is  a  quantity  which  is  negligible  in  planar 
approximation.  (Moritz,  1969). 


In  most  cases,  the  cap  size  will  be  sufficiently  small  to  replace  (3.6)  by  a 
plane  integral. 


In  that  case  we  proceed  in  two  steps .  We  assume  that  —  gp(5 )  has  to  be 

R 

evaluated.  In  doing  this  we  replace  the  integration  over  the  unit  sphere  by  one  over 
a  sphere  with  radius  P,.  (3.6)  becomes  then  after  neglecting  the  last  integral 


RgC(?)  l  IrV-’r^P  R  dT(r,) 


=  I  {  — — -  ■-  -  -  - - -  }  -5-  Lapu<p(§,7j)R2dr(77) 

c  l  IR5  -  Rrjl  R  je(cos  0q)  J  R1  ™ 

cos  & 

-  - <p(5,  r\)  R  d  aC(rj) 

(Rf(cos  4>o)r*c 


(3.7) 


Now  we  make  the  transition  to  the_plane  replacing  RgJjy  x  =  (-xj,X2)  » 

T  1 

R77  by  y  =  (yp  y2>  and  denoting  by  p  the  radius  of  the  cap.  The  expression  ^ 

Lap qViZ'W)  goes  over  into  the  ordinary  plane  Laplacian  lapy<p(x,y)  =  <Py^y  +  ^2^2’ 


Hence: 


11 


lyl 


=  1  {ttV  '  p)  “vy  v(x-y)  dy 

lyl 

2n 

-  -  f  <p(x,x  +p  /  008  a  |  )  da  (3.8) 

p  ^  l  sin  a  J 

a=0 


It  will  be  clear  from  later  sections ,  that  we  shall  have  a  polynomial  re¬ 
presentation  for  the  function  <p(x,y).  However,  this  polynomial  representation  will 
vary  from  mesh  to  mesh  of  a  rectangular  grid.  It  may  therefore  be  desirable  to  have 
a  cap  boundary  which  is  not  circular  but  rectangular.  Call  this  rectangular  area  Q. 
Then  (3.8)  becomes: 


; 


lx  -  yl 


kPy<p(x,y)  dy 


+  <j>  <p(x,y  )  (grady  --  -*-j- ,  v(y))ddQ(y) 

?)Q  X’y 


P  — - —  (grady<p(x,y),  v(y))ddQ(y)  (3.9) 
SQ  lx-yl 


These  formulas  could  also  be  verified  by  Green’s  second  formula  for  the  plane  consider 


ing  that 


laPv 


lx-yl 


|x  -  yf 


(3.10) 


v  is  the  outer  normal  to  the  boundary  of  Q  and  grad  denotes  the  ordinary  grad 
icnt  in  the  plane.  Thus  the  second  equal  sign  in  (3.8)  and  (3.9)  is  exact.  Only  the 
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equality  with  the  integral  over  the  portions  of  the  sphere  is  approximate. 


3.3  Reg  uiization  of  Vening -Meinesz '  formula 
Vening  Meinesz'  formula 

v<5)  =  -  fGrad  St(M?)Ag(n)dIXr?)  (3.11) 

constitutes  also  a  singular  integral  since  it  does  not  exist  after  replacing  the  integrand 
by  its  absolute  value.  Regularization  may  take  place  in  a  similar  way.  Again,  only 
the  contribution  v^(  5)  of  a  circular  cap  C  centered  at  §  and  having  spherical  rad¬ 
ius  il>0  is  considered.  Farther  outside  the  original  form  (3. 11)  is  quite  suitable. 

The  regularization  procedure  is  outlined  in  Appendix  D.  It  results  in: 

vc<?)  =  jffG  Ja(?’T7)Lap^r?)dr<T?)  ‘ 

C 

-  r~  §  o(S,n)  (Grad  Ag  (77) ,  v(n ) )  d a  C(n ) 
zIT  & 

ac 

3  $ 

+  27c  !  - ~~7T2'  cr(§,  T7)  >  c! r  (77  )  - 

C 

-  |  Grad?R(5-rj)  Ag(r7)dT(77) 

C 

Thereby  we  have  split: 

st(  n )  =  -r—T  +  R  (?•»?)  (3.13) 

£(C-r)) 

The  unit  vector  a(£,  77)  is  tangential  to  T  in  the  point  §  and  lies  in  the  same 
plane  as  *  and  77.  R (5-77)  causes,  even  after  differentiation,  no  singularity.  The 
singularity  in  Vening  Meinesz’  formula  comes  in  through  the  term  2/£ ( § *77)  in  (3.13) 
and  has  been  removed  in  the  above  formula.  In  Appendix  D  there  is  also  a  listing  of 
the  above  formula  in  the  usual  notation  which  is  based  on  a  localized  coordinate  system 


+ 


(3.12) 


13 


(in  which  a  (£  •  r))  carries  over  into  (  gin  ^ ) .  Equation  (D.  9a). 

Again,  formula  (3. 12)  would  become  simpler  if  C  would  be  replaced  by  T.  The 
line  integrals  would  vanish.  The  numerical  Usefulness,  however,  would  decrease  for 
reasons  similar  to  those  given  in  section  3. 2. 


3.4.  Some  further  manipulations 

There  is  still  a  singular  integral  in  the  deflection  formula  (1.2).  It  is 


q(£) 


4ff  R2  G 


P  _  2  1 

1  ( h  (r> )  -  h  (£  )  Grad  '  3  " 


£  l  (t-V) 


Ag(T?)  df  (Ti) 


This  term  is  combined  with  a  certain  other  contribution  toward  (1. 2)  and  regularized 
as  outlined  in  Appendix  E.  A  summary  of  all  regularizations  is  given  in  the  following 
subsection. 
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3.5  Summary  of  regularization 

3.5.1  The  Quantity  G^(5). 

For  the  area  outside  the  cap  C  use  the  original  formula  (1.3).  For  the  cap- 
interior  use: 

G.  (|)  =  -JL_  J  -  - I -  }Lapr,[(h(rj)-h(5))2ig(r?)]dr(77) 

LC  2irR  c  1  Hcos  J  I 


§  (hfo)-h(g)  )Ag(n)da  C(T7)  (3.18) 

£2(cos  0Q)  ac 

This  formula  follows  from  (3.6)  after  neglection  of  the  last  integral  (planar  approximation). 

3.5.2  The  quantity  G2(f) » 

Split  __ 


G2(§)  =  G2(5)  +  G2(§) 

(3.19) 

C2<5) 

-  G1(i))dr(ij) 

/.IT  n''  /  ~  \ 

r  *  (£•??) 

(3.20) 

S2<5» 

=  ~2  Ag(?)  |  Grad  h(?)|  2 

(3.21) 

For  G2r  use  (3. 18)  with  Ag(rj)  replaced  by  G  1(77).  The  evaluation  of  G^^)  is 
immediate . 


3.5.3  Formula  (1.1)  For  The  Height  Anomaly  a(g). 
No  further  regularization  is  necessary. 
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3.  5. 4.  The  Deflection  -  Formula. 


For  the  area  outisde  the  cap  C  the  original  formula  (1.2)  will  be  used. 

For  the  cap  interior  the  following  will  be  used  (cf.  (3. 13) ): 

Vc<?)  =  ~dc  $  CT(5*’7>  LapCAfKrj)  +  G1(tj)  +  G2(t))]  dT(tj)  - 

c 

■  £  c(S.  r?)  (GradLdgftj) +G,(r)) +G2(t))]  ,  u(t)))  d5C(?))  + 

2  TT  G 

3C 

.  3$ 

,  „  1  -  COS  Tf 

+  — L.  j*  - a(?,T))[Ag(T7) +0,(77) +G2(n)l  dr  <1|>  - 

2WG  sinzib 

C 

-  J~  j*  Grad?R(5.r?)[^g(T7)+G1(t7)+G2(r))+B2(r7)]  dr(7))  + 


4r8*G  J  ®<M>  { I«P,[(l>('!)-h(j))2Ag(7)] 

-  2|Grad  h(^) J  2Agty) j  d  TC^) 


■5  v 

r  1  -  cos-'  -w 

- 5“  J  - : — j - 

2  tB  5  ^ 


0-(j»7)  (-Grad  h(i|)  |  2  Ag(r| )  d  T  ( 17) 


— L—  f  — V  crU, 7)  (  ffradJ(h(7)-h(f))2Ag(7), 
V  ir  P  G  Bin  y  ‘  J 

»  V*  (7)  )  d&IO?) 


if  o  cos  V" 

- f  [(h(7)-h(|))2  AgC^)]  -  -t  cr(f  ,7)  d  2> C (7) 

2  ir  r.  G  ££  sin  iy 


-1  IT  K' 


J  0(>))-h(!))2Ag(7)] 


3(1  -  cos-'  -p) 

- — j - - 

sin  y 


■  C3'(I»7)  dP(7)  (3.22) 
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4.  ANALYTIC  REPRESENTATION  OF  FUNCTIONS 

It  has  been  sated  repeatedly  that  our  regularization  procedures  are  based 
on  the  availability  of  an  analytical  representation  of  the  various  involved  functions 
like  terrain  heights,  gravity  anomalies,  and  other  functions  derived  from  them.  Data 
libraries  contain  information  about  these  functions  either  in  discrete  form  or  in  block 
average  form.  Neither  of  these  representations  is  suited  for  our  purposes.  Block  aver¬ 
ages  produce  step  functions  with  constant  values  within  each  block  and  jumps  at  the 
block  boundary.  Such  a  function  is  not  analytical  nor  can  it  be  reconstructed  from  its 
first  and  second  derivatives  which  are  zero  almost  everywhere.  For  our  regulari¬ 
zation  procedure  it  is,  however,  essential  that  the  functions  can  be  reconstructed 
from  these  derivatives  (plus  some  additional  information  like  the  function  values  on 
a  boundary  line). 

If  gravity  anomalies  and  terrain  heights  are  given  at  discrete  values,  an  inter¬ 
polation  procedure  has  to  be  used.  For  an  irregular  pattern  of  sample  points  in  two 
dimensions  this  can  be  a  quite  laborous  task.  Usually  a  linear  procedure  is  used 
leading  to  an  expression 

N 

P(x,y)  =  £  ziPi(x,y)  (4.1) 

i  =  1 

z^  are  the  values  of  the  functions  for  the  arguments  Xj,  yj,  p(x,y)  is  the  interpolat¬ 
ing  function.  The  functions  P|(x,y)  depend  on  only  the  arguments  x^,  y^  and, 
of  course,  on  the  chosen  interpolation  procedure.  They  do  not  depend  on  the  z^. 

Examples  are:  polynomial  interpolation?trigonometric  interpolation,  linear 
prediction.  If  the  function  Pj(x,y)  is  differentiable,  then  p(x,y)  is  differentiable. 
This  generally  solves  our  problem.  The  question  is  how  to  choose  a  computationally 
feasible  interpolation  procedure. 

For  an  irregular  pattern  of  points  there  is  hardly  a  chance  to  get  simple  ex¬ 
pressions  for  p(x,y).  Therefore  a  rectangular  grid  of  points  (x^,yj)  points 'must 
he  achieved.  This  can  be  done  for  example  by  a  prediction  procedure.  The  interpola¬ 
tion  within  the  rectangular  grid  is  according  to  the  author ’s  opinion,  in  the  best  way 
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accomplished  by  spline  interpolation.  Spline  interpolation  yields  a  polynomial  expres¬ 
sion  for  the  function  within  each  mesh  xi  ^  x  £x^+^  yj  <.y  of  the  grid.  The 

polynomial  expression  varies  from  mesh  to  mesh.  However,  continuity  and  continuous 
differentiability  up  to  any  desired  order  can  be  maintained  at  the  mesh  boundaries. 
Advantages  of  spline  interpolation  over  other  methods  are  outlined  in  Appendix  A. 
Loosely  speaking,  spline  interpolation  keeps  the  derivatives  of  the  interpolating  func¬ 
tion  as  small  as  possible,  thus  avoiding  fluctuations  and  oscillations  as  they  are  fre¬ 
quently  encountered  in  other  methods.  Another  advantage  is  that  the  polynomial  repre¬ 
sentation  within  a  certain  mesh  depends  to  a  large  extent  only  on  the  function  values 
y.-j  in  the  near  neighborhood  of  the  mesh.  The  dependency  on  values  farther  outside 
is  negligibly  small.  This  is  also  not  the  case  with  other  methods  where  a  change  of 
one  value  z^  may  cause  changes  in  the  interpolating  function  over  a  wide  area. 

in  Appendix  A  computational  procedures  are  outlined  to  obtain  the  interpolat¬ 
ing  bicubic  spline  functions.  Asymptotic  formulas  have  been  used  which  are  simple 
and  quite  appropriate  for  our  purposes. 

Spline  interpolation  can  also  be  applied  toward  the  block  averages .  It  is 
only  necessary  to  view  the  block  averages  in  a  different  way.  They  should  not  be 
viewed  as  step  functions.  Instead  a  block  average  shall  be  representative  only  for 
its  midpoint.  For  any  other  point  the  appropriate  value  would  be  an  average  over  a 
block  centered  at  that  point.  In  other  words  we  imagine  a  smoothed  version  of  our 
function  which  is  obtained  by  moving  averages  over  blocks.  Unfortunately  we  have 
values  of  this  smoothed  version  only  at  discrete  values,  i.e.  the  midpoints  of  the 
original  blocks.  For  the  other  points,  an  interpolation  procedure  has  to  be  used. 

This  can  certainly  be  spline  interpolation  since  the  midpoints  already  form  a  rectang¬ 
ular  grid. 

Speaking  of  a  rectangular  grid  does  not  mean  that  the  lines  must  be  equidistant 
on  the  unit  sphere  f.  The  grid  lines  can  certainly  be  represented  by  lines  of  constant 
latitude  and  longitude.  In  evaluating  the  various  differential  operators  like  Grad  and 
Lap  the  appropriate  formulas  have  to  be  used. 


Spline  interpolation  will,  of  course,  only  be  used  for  the  small  cap  areas  in 
which  regularization  takes  place.  For  the  outer  zones  the  original  averages  or 
even  cruder  representations  will  be  sufficient.  Here  one  benefit  of  spline  interpolation 
becomes  effective.  The  polynomial  representations  for  a  mesh  is  practically  depend¬ 
ent  only  on  die  function  values  in  a  rather  small  neighborhood  of  the  mesh.  This 
avoids  the  necessity  of  performing  an  a  priori  interpolation  for  larger  areas  where 
continuity  problems  arise  at  the  connecting  boundaries. 
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TRUNCATION 


3  . 

The  procedure  for  evaluating  the  various  integrals  will  be  the  following.  For 

AS 

a  very  small  cap  with  half  opening  angle  0^  the  regularized  formulas  will  be  used 
in  cast-  the  integral  is  singular.  The  radius  of  this  cap  will  be  a  few  kilometers  or 
even  less.  This  corresponds  to  a  0^  of  about  1‘.  Outside  this  cap  essentially  the 
original  form  of  the  integrals  will  be  used.  It  is,  however,  certainly  not  necessary  to 
carry  out  the  integration  with  the  same  accuracy  all  over  the  sphere.  Some  integrals 
can  be  completely  truncated  at  a  certain  distance.  Some  have  to  be  extended  over 
all  of  1'.  However,  going  farther  out  more  smoothed  versions  of  the  functions  h, 

Ag,  . . .  can  be  used.  This  reduces  the  computation  time. 

There  exists  a  theory  on  truncation  which  goes  back  to  Molodensky.  See 
Molodensky  et  al.  (1962).  De  Witte  and  others  have  done  further  work.  In  Appendix  B 
I  have  reviewed  and,  as  I  think,  somewhat  extended  this  theory. 

New  insight  has  been  gained  into  the  nature  of  the  truncation  error.  Take  for 
example  Stokes'  formula.  De  Witte's  work  suggests  that  truncation  is  most  favorable 
at  the  zeros  of  Stokes'  function.  However  a  zero  can  be  placed  anywhere  by  simply 
adding  a  constant  to  Stokes'  function.  Doing  this  and  truncating  at  such  an  enforced 
zero  produces  a  genuine  truncation  error  which  is  comparable  in  size  to  that  after 
truncation  at  a  zero  of  the  unmodified  Stokes’  function.  The  only  difference  is  that 
the  new  truncation  error  is  superimposed  by  a  constant  which  is  the  same  for  all  points 
of  the  sphere.  Such  a  constant  does  not  matter  in  Stokes'  formula.  The  geoid  has  to 
be  scaled  later  on  anyway. 

We  have  also  investigated  the  effect  of  replacing  the  abrupt  truncation  at  a 
certain  angle  by  a  smoother  procedure.  It  consists  of  leding  the  kernel  taper 
off  to  zero  over  a  certain  interval  0j  <;  0  ^  0q  •  This  brings  along  a  very  beneficial 
effect  in  dampening  the  higher  harmonics  of  the  truncation  error.  The  situation  is 
similar  to  filter  design  theory.  Filter  functions  with  discontinuities  tend  to  cause 
errors  wilh  considerable  high  frequent  portions.  Smoothing  out  the  discontinuities  re¬ 
moves  these  high  frequent  errors  to  some  extent. 


20 


In  the  following  we  discuss  truncation  procedures  for  the  various  integrals. 
Alternatives  to  the  proposed  strategies  can  be  obtained  from  Appendix  B. 

5.1  The  integrals  of  the  correctional  type 

Gj  is  a  correction  of  Ag  as  it  is  seen  from  (1. 1)  and  (1.2).  Therefore  a 
contribution  toward  Gj  can  be  neglected  if  it  amounts  to  less  than  the  measuring  ac¬ 
curacy  of  Ag.  This  accuracy  may  be  assumed  as  one  part  in  1000.  Due  to  this  and 
also  due  to  the  fact  that  the  function  l/£  ^  ( 5*  rj )  tapers  off  quickly  truncation  can 
take  place  at  a  short  distance.  Bursa  (1965)  proposes  truncation  beyond  a  distance  of 
about  80  km.  This  corresponds  to  ipg&O.  7°.  A  global  statement  is  difficult  to  make 
but  some  rough  estimates  based  on  the  truncation  theory  of  Appendix  B  indicate  that 
even  in  the  case  of  our  second  order  approach  not  much  more  is  needed.  i/)q«=>  2 °, 
should  be  sufficient  in  all  cases. 

The  same  can  be  said  about  all  other  integrals  of  this  correctional  type. 

These  are:  The  integral  in  G2>  the  second  integral  in  either  of  (1. 1)  and  (1.2) . 


5 . 2  Stokes '  integral 

Whereas  nothing  new  nas  been  said  in  section  5.1,  it  is  hoped  that  the  following 
discussions  on  Stokes'  and  Vening  Meinesz  integral  bear  some  interest. 

We  propose  the  following.  Split  off  the  spherical  harmonic  components  up  to 
and  including  degree  12  from  Stokes'  kernel: 

12 

St(l'Tj)  =  £  pn(?*r?)  +  St  12(§* r?)  (5.1) 

n  =  2  n "  1  n  u 

The  contribution  of  the  first  part  will  not  be  discussed  further.  It  could  be  taken  from 
satellite  solutions.  The  residual  kernel  St  ^( cos  $)  has  at  least  13  zeros  in 
0  ::  1 1)  180.  (  Meissl  (1971)  )  .  The  zeros  are  approximately  given  by 

ll)  -  5.02°,  17. 99°,  31.58°,  45.29°,  59.06°,  72.85°,  86.65°,  100.5°, 
114.3°,  128.1°,  141.9°,  155.7°,  169.4°. 
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These  zeros  are  good  locations  for  truncation,  or  better  for  a  transition  to  a  heavier 
smoothed  version  of  Ag  in  the  residual  Stokes’  formula: 


n12(5)  -  -5-  J  st12<s-t))  4g(,)dr<,> 

4tro 


(5.2) 


In  (1. 1)  Ag  is  replaced  by  ig  +Gj  +G2  •  However  for  the  sake  of  simplicity  we  dis¬ 
cuss  the  (residual)  Stokes  formula  in  the  form  (5.2). 

Now  we  integrate  only  up  to  \1)q  =  5.02  with  full  accuracy.  Beyond  $  =  5.02 
we  replace  Ag(§)  by  a  smoothed  version  Ag(  ^(4)  which  is  obtained  by  averaging 
Ag(§)  over  a  circular  disk  of  half  opening  angle  0^  centered  at  5.  The  spherical 
harmonic  expansion  of  the  resulting  error  AN^(?)  in  N12(f)  is  given  by  (B.82) 

and  (B. 45)  (0) 

foi  r  Qn  (0) 

7T  -T-d-rf  <5*3> 


(0) 

(AN,) 

12  nm 


^  are  the  Molodensky  coefficients  referring  to  the  residual  Stokes '  kernel 

St12(  §'t?)  and  a  truncation  angle  equal  to  .  The  Q„°  ^  can  be  computed  from 
(B.46)if  St( 77  )  is  replaced  by  Sti2(§*T))  and  the  summation  over  m  is  extended 
to  N  =  12.  The  ^  are  the  eigen  values  of  the  smoothing  operator  over  the  Oq  - 
disks.  The  p^  ^  are  given  by  (B.83)  or,  for  not  too  large  n,  approximately  ly 
(B.86). 


Let  us  use  the  smoothed  version  Ag(  * )  only  up  to  =  17.99.  CXitside  the 
ib  j  cap  a  still  heavier  smoothed  version  Ag^  \  will  be  employed.  Ag^^  is 
ohained  by  averaging  Ag  over  disks  of  half  opening  angle  dtp  >  oi^.  The  total 
error  would  then  be  given  by 


(AN*;*)) 


nm 


(1) 

nm 


>Agnm 


(5.4) 


Molodensky  coefficients  for  truncating  St ^2( c<>s  0)  at  0=$^. 
eigen  values  of  the  averaging  operator  over  the  0^  -  disks. 

This  procedure  is  immediately  generalized  to  a  succession  of  zones.  If 
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is  the  last  truncation  angle  then  the  associated  error  in  Nj^  Is  AN  12*  The  following 


recursion  holds: 


nV1  / 


The  general  form  of  the  error  is 

/AM^) 


(AN*  *2  )  =  ^  ^  Ag 

12  inn  O  n  nm 


2  o  (i  ) 

If  a  ( Ag)  are  the  degree  variances  of  Ag  then  the  mean  square  error  a  (AN  J2  ) 

is  given  by 


a2(AN(i2))  =  £  (q(*  V  a2(Ag) 

12  G2  n>12  n  n 


(5.7) 


See  for  example  Meissl  (1971)  chapter  4,  esp.  equ.  (4. 11)  there., 

We  adopt  the  following  degree  variances  for  Ag  (they  are  taken  from 
Gaposchkin-Lambeck  (1970).  ). 

2  2 

n  an(Ag )  in  mg 


Occasionally  we  shall  need  higher  degree  variances  too.  Since  we  need  them  only  for 
rough  estimates  we  adopt  the  following  model  for  the  <j2(Ag),  n  >  12: 


(5.8) 


CTn(Ag)  = 


3.10 


11 


(n  +  445 f 


Summing  all  degree  variances  should  give  a  value  close  to  1201  mg  .  (Kaula  (1959). 
We  have  from  the  above  table 

12  2 

E  o*(Ag)  =  159  (5.9) 

n  =  2 


and  from  (5.8) 


E 

n>  12 


a2(Ag) 


1042 


(5. 10) 


o 

which  gives  together  1201.  Besides  that  we  have  for  n  =  13, a  value  a  (Ag)  s»  7  mg 

13 

which  fits  into  the  picture  of  the  above  table. 


The  variance  of  the  geoidal  undulations  is  then 


2 

a 


E  — 9  CT2(4g)  =  (30.2m)2 

n  =2  (n-D2  n 


(5.11) 


This  is  the  total  variance.  The  variance  of  the  residual  part  is 

2 

a2(N,J  =  E  — L_„  a2(Ag)  =  (6.3m)2  (5.12) 

12  G2  „>12  (n-1)2  n 

(i) 

In  Table  5-1  we  have  listed  q  quantities  under  various  assumptions  con- 

n 

ceming  the  underlying  o^’s. 

Take  for  example  an  exact  representation  of  Ag  within  0 q  =  5.02  and  moving 

averages  over  circular  disks  with  half  opening  angle  a «  =0.56  outside.  This  corres- 

(0)  v 

l>onds  roughly  to  1°  x  1°  blocks  outside.  The  qn  are  seen  not  to  exceed  0.0001 
for  n  >  12.  This  enables  us  to  compute  an  error  bound  which  is  independent  of  the 
hypothesis  (5.8).  From  (5. 7)  follows : 

a2(AN(°2))=  A!  E  (q(0))2a2(^)  s  A!  0.0001  E  a2(Ag) 

12  G2  n>12  n  n  G2  n>l2  " 
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From  (5. 10)  follows  now 


<?(AN(0))  £  JL  0.0001  Vl042'  **  0.02 
12  G 


(5.13) 


Under  the  hypothesis  (5.8)  even  a  smaller  error  can  be  estimated:  a(AN 0.006  m. 

Truncating  further  ac  =  17.99  and  using  moving  averages  based  on  =  1.4 

outsioe  (roughly  2.5°  x  2.5°  -  blocks  )  produces  an  error  not  larger  than  .11  m.  This 

again  independent  of  hypothesis  (5.8).  This  error  would  be  obtained  if  all  the  power 

o 

of  Ag  above  degree  12  was  concentrated  at  n  =  13:  a  (Ag)  -  1042.  This  is  certainly 

13 

unrealistic.  Under  hypothesis  (5.8)  the  estimate 


a(AN^)  0.02  m 


(5.14) 


holds . 


Truncating  further  at  $  =  31.58  and  using  =  2.8  moving  averages  out¬ 
side  (roughly  5°  x  5°  -  block  avei-ages)  gives  a  total  error 


(2) 

o(AN  l2  )  «  C. 33  m 


(5.15) 


Hypothesis  (5.8)  has  been  used  thereby. 

5.3  Vening  Meinesz's  integral 

Again  we  split  according  to  (5. 1)  and  assume  that  the  contribution  due  to  the 
harmonics  up  to  and  including  n  =  12  can  be  taken  care  of  otherwise.  If  truncation 
takes  place  at  the  zeros  of  St  (cos  $  )  ,  then  the  theory  of  Appendix  B,  section  B.6, 
tells  us  that  the  variance  of  the  resulting  error  in 

v12(?)=  — —  J*  Grad  Stl2(5'T7)  Ag(7?)dr(tj)  (5.16) 

r 

is  given  by 

a2(Av(,V)=  -V  £  n(n+l)  (qj°)2  a2(Ag )  (5.17) 

G  n>  12  n 

(O 

See  erju.  (15.  07a).  The  qn  are  those  of  the  previous  section. 


23 


Table  5-1 


The  below  are 

n  * 
given  in  unite  of  10“°. 

If  the  n-oolumn  shows 

an  interval  rather  than 

a  number, then  the 

extremum  values  within 

this  interval  are  listed. 


n 

5(°) 

q(1) 

*n 

*n 

13 

91 

523 

1905 

14 

90 

446 

1314 

15 

86 

361 

739 

16 

85 

274 

251 

17 

82 

189 

-100 

18 

79 

110 

-290 

19 

75 

42 

-328 

20 

70 

-13 

-250 

21 

66 

-52 

-108 

22 

60 

-75 

42 

23 

55 

-82 

151 

24 

49 

-75 

193 

25 

44 

-57 

165 

26 

38 

-30 

85 

27 

32 

1 

-15 

28 

26 

32 

-98 

29 

20 

60 

-138 

30 

14 

82 

-120 

31-  35 

-13 

100 

226 

c 

'j- 

i 

-33 

-91 

-218 

41-  50 

-44 

-118 

-254 

51"  73 

-42 

74 

-142 

76-100 

34 

64 

98 

101-150 

-28 

-55 

-64 

151-200 

24 

44 

44 

201-300 

18 

20 

20 

3°l-300 

12 

11 

11 
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The  overall  variance  of  the  deflection  of  the  vertical  is  given  by 


a2(v)  =  -L-  z  q2(Ag)  (5.18) 

G 1  n^2  (n-1)2  n 

under  the  hypothesis  (5.8)  this  leads  to 

a(v )  =  7.9"  (5.19) 

The  residual  deflections  v^  have  standard  deviation  obtainable  by  summing  in 
(5. 18)  only  over  n  >  12.  This  gives 

a(v12)  =  6.9"  (5.20) 


The  truncation  procedures  outlined  in  the  previous  section  produce  errors  which  can 
be  estimated  in  the  following  way. 


*0 

-  5.02. 

oc  o  0  •  56 

o(Av<12))  *  0.03" 

(without  hypothesis  (5.8)) 

*1 

=  17.99, 

0^  =  1.4 

a*  0.013" 
a(Av<l2))  *  0.02  " 

(under  hypothesis  (5.8)  ) 

(under  hypothesis  (5.8)  ) 

*2 

=  31.58, 

«2  =2.8 

a(Av(2>«  0.03" 

12 

(under  hypo  diesis  (5.8)  ) 

03 

=  45.29, 

a3  =5.6 

^(Av^s*  0.07" 

(under  hypothesis  (5.8)  ) 
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6. 


A  COMPUTATIONAL  TEST 


I  am  fully  aware  that  any  proposed  numerical  approach  toward  a  mathematical 
problem  should  be  thoroughly  checked  by  test  calculations.  Unfortunately,  time  did 
not  permit  to  do  tiiis  on  a  large  scale  for  the  present  study.  In  order  to  give  some 
idea  how  the  proposed  methods  can  be  converted  into  computer  programs  I  shall  deal 
in  the  following  with  bicubic  spline  interpolation  and  with  the  -integral.  We  assume 
a  plane  grid  which  is  quadratic  and  where  the  spacing  between  the  lines  is  unity.  We 
do  not  bother  about  a  realistic  scaling  of  the  model.  Let  the  grid  consist  of  points . 
Assume  the  functions  hjj  and  Ag^  be  defined  for  the  grid  points  by 


N 

hij  =  E 


k  =  l  i  +2j  +  3k 


sin(i~.)  sin(  lX?L) 


N 


Ag  ■  =  £  - - 

U  k  =  1  i  -t-  2.5  j  +  4k 


sin(  sin(JJ^L) 

N  N 


(6.1) 

(6.2) 


After  this  data  generation  we  proceed  to  the  bicubic  spline  interpolation.  Let  us  out¬ 
line  this  for  h.j .  The  procedure  for  Ag^  is  completely  analogous.  Compute  ac¬ 
cording  (A.  7) 

a  =  -  (2  -  \(3)  -.268 


Now  according  (A.  16)  to  (A.  19)  compute  with  0  =  6 
6 


m. .  = 
ij 

-  3 

£ 

k-1 

i  =  7, ....  N-6 

j  —  1, . .  • ,  N 

(6.3) 

6 

nij 

-  3 

£ 

^^i,  j+£  "  hi,  j -P 

i  =  1, ...,  N 

(6.4) 

j  =  7,...,  N-6 

6 

I1  = 

-  3 

£ 

c^(m.  .  ,  „  -  m.  .  .) 

i,j+£  i,j-J 

i  =  7.....  N-6 

ij 

t  =  1 

II 

? 

a- 

(6.5) 

A  check  could  be  made  by  using 
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MW***' 


p.  =  -  3 
ij 


S  o  ( rij  +jcj  “  nj  .  ^  j) 

k  =  1 


i  =  7, ....  N-6 
j  —  7,  . ..,  N-6 


(6.6) 


if  the  agreement  with  (6.5)  should  not  he  satisfactory,  then  a  /3  larger  than  6  has 
to  be  chosen.  In  our  test  /3  =  6  was  sufficient.  Now  h^  can  be  interpolated  for  any 

one  of  the  meshes  i  £  x  s  i  +  1,  j  ^  y  £  j  +  1  within  the  subgrid  i  =7 . .  N-7, 

j=7,  _ _  N-7.  The  formulas  are 


MX’y)  ‘  kjl  =0 


i  x  £  i  +  1 

j  s  y  ^  j+1 


(6.7) 


,(i»j  ) 


The  a  j,  g  ’  are  computed  with  the  aid  of  (A.  13)  to  (A.  15).  One  can  pre-compute  the 

u^’-i )  for  as  many  meshes  as  storage  allows.  They  can  always  be  computed  from  the 
k,  C 

h:;,  niM,  nfiJ  p:  ,  which  do  not  require  that  much  storage  (one-fourth). 

AJ  AJ  AJ  A>J 

Turning  now  toward  the  G  j-  integral  take  a  plane  approximation  to  formulas 
(1.3)  and  (3.  18): 


G1(x,y)=— L 

2  7T 


S 


d5r  d?  (6.8) 


(x-x)2  +  (y-y)2  ^  R2 


[(x-x)  +  (y-y)  ] 


1 


Glc(x,y)=^  2  1  2  •  •  1 

(x-x)  +  (y-y)  *  p 


•  iap— -  C(h(x,y  )-h(x,y  )  Ag(x,y  )]dx  dy  - 
» y 


i  i 
2 n  p"2 


Sj  [h(x,y  )-h(x,y  )]  Ag(x ,y  )  d  A  C(x,y  ) 

AC 


Converting  everything  to  polar  coordinates  centered  at  x,y  gives: 


(6.9) 
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(6.10) 


R  2  it 

Gl(x,y)=^  J*  J  Oh(x +r  cosa,  y +r  sina)  -  h(x,y)] 

r  =0  a=0 

j 

•Ag(x+r  cosa,  y  +  r  sina)  da 

r  *- 

p  2tr 

Glc<X,y>=  2?  J-  -T  [1'p]liIp3!,y  £W5T.m<*.y>^<S,y%x+rcosa 

r  =0  a=0  y=y+rsina 

.  dr  da 

2  it 

1  1  ,  r  -i 

■  —  Lh(x  +  r  cosa,  y  +  r  sina)  -  h(x,y )  J 

a=0 

Ag(x+r  cosa,  y+r  sina)  da  (6.11) 

I  have  used  the  following  procedure  for  the  evaluation  of  the  integrals .  In  (6. 10)  the 
integration  with  respect  to  a  has  been  carried  out  first  using  the  trapezoidal  rule. 

(The  trapezoidal  rule  is  advantageous  for  periodic  functions.)  Integration  takes  place 
for  a  certain  set  of  discrete  r-values .  Afterwards  the  trapezoidal  rule  has  been 
applied  to  integrate  over  r.  Values  for  h  and  Ag  have  been  computed  from  spline 
formulas  throughout.  (This  has  only  been  done  for  checking  purposes.  Later  on  we 
will  use  a  less  complicated  procedure  for  the  more  distant  areas  )  The  first  integral 
in  (6. 11)  has  been  treated  similarly.  The  expression  involving  the  Laplacian  has  been 
expanded  as: 

Jajv-  - [(h(x,y)-h(x,y) )  Ag(x,y)]  =  [h(x,y  )-h(x,y)]  lap  Ag(x,y) 

+  Ag(xVy)  laph(x,y)  +  2  (grad  h(x,y),  gradAg(x*,y)  (6.12) 

The  various  derivatives  have  been  computed  from  the  spline  interpolation  formulas. 

(Note  that  the  bicubic  spline  has  continuous  second  derivatives.)  The  second  integral 
in  (6.  11)  has  also  been  approximated  with  the  aid  of  the  trapezoidal  rule. 

For  N=21  we  took  x=ll,  y=ll,  p=.8,  R=3.9.  The  result  was 
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g2(  11, 11)  =  - 


—  0. 1525 
2  IT 


For  checking  purposes  we  also  took  p=.  5  .  The  result  was  then 


Gi(  11. 11)  =  -  ~  0.1510 

1  2  it 


The  various  stepsizes  for  the  numerical  integration  were 


,  .  P  .  2  IT 

for  r  s  p  :  Ar  =  —  ,  Aa  =  - 

8  4  r/Ar 


k-1 


for  r  s  p  :  Ar  =  const.  1.3  ,  k  =  1, . . .,  8,  const 


Aa  =  2  tr/32 
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APPENDIX  A:  Some  Aspects  of  Spline  Interpolation 

It  is  not  our  intention  to  give  a  systematic  accountof  spline  interpolation.  Books 
like  Ahlberg  et  al.  (1967)  may  be  consulted  for  a  thorough  discussion.  The  idea  be¬ 
hind  spline  interpolation  and  some  of  its  advantages  shall  be  outlined.  Formulas  and 
methods  used  in  the  main  portion  of  this  report  shall  be  documented. 

A.l  The  Cubic  Spline 

Let  us  start  with  a  sequence  of  n  equidistant  abscissa's  xi#  i=l,  ...,n  on 
the  x-axis  of  a  two-dimensional  coordinate  svstem.  Let  the  distance  x.  ,  .-x,  be 
denoted  by  h.  Let  function  values  y^  be  prescribed  for  each  Xj  .  The  old  problem 
of  interpolating  these  discrete  function  values  can  be  solved  in  many  ways.  One  way 
is  interpolation  by  fitting  a  polynomial  of  degree  n-1  through  the  n  points.  This 
can  be  done  e.g.  using  a  polynomial  p^(x)  of  degree  n-1  having  the  properties 
p.(Xj)  =  1  for  i=j  and  p^(xj )  =  0  otherwise.  Analytical  expression  for  these 
polynomials  are  found  under  "Lagrange  interpolation"  in  any  relevant  textbook.  The 
interpolating  polynonv  .1  is  then  given  by 

n 

p(x)  =  E  (A*1) 

i  =  1  1  1 


An  advantage  is  the  analyticity  of  the  resulting  expression  which  is  infinitely 
differentiable.  A  disadvantage  is  the  following.  A  charge  (an  error)  in  one  of  the  y’s 
causes  considerable  changes  (errors)  over  most  of  the  whole  range.  This  is  immed¬ 
iately  seen  from  equation  (A.l).  A  change  Ay,  in  y.  causes  a  change  in  p(x  ) 

Ao  *o 

which  is  given  by 

Ap(x)  =  Ay,  p.  (x)  (A. 2) 

o  *o 

This  may  not  be  too  bad  as  far  as  only  function  values  are  concerned.  The 
changes  may  be  much  larger  for  the  derivatives.  An  extremely  opposite  approach 
would  be  piecewise  linear  interpolation.  The  interpolated  function  is  then  between 
x,  and  x  given  by  the  straight  line. 

1  i  +  1 


X  -  Xj 


J2(x) 


yi 


X  -  Xj  +  1 
Xj  -  Xi+1 


yi+i 


xi+r 


(A.  3) 


A  change  in  y  causes  now  only  changes  in  Xj-  i  £  x  £  xi  +  j  and  no  further 
reverberations.  A  disadvantage  is  evident  from  the  fact  that  the  function  is  not  con¬ 
tinuously  differentiable  over  the  entire  range. 

Spline  interpolation  tries  to  compromise  between  these  two  extremes.  The 
interpolating  function  is  assumed  to  be  a  polynomial  of  odd  degree  k  in  each  of  the 
intervals  Xj  <  x  <  x^  +  j  .  However,  these  polynomials  vary  from  interval  to  inter¬ 
val.  The  continuity  of  the  function  as  well  as  of  its  first  k-1  derivatives  is  required 
throughout,  i.e.  also  for  the  abscisses  x^,  i=l, .  ..,n.  Certain  boundary  conditions 
for  xj,  xfl  may  be  prescribed  in  order  to  make  the  problem  unique.  This  interpolat¬ 
ing  function  s(x)  is  then  (k-1)  times  continuously  differentiable  and,  as  we  shall 
see  later,  a  change  in  one  of  the  y.  causes  changes  in  s(x)  which  are  heavily  damp¬ 
ened  for  larger  distances  x  -  v^. 

We  discuss  here  only  the  case  k=3  which  is  known  under  the  label  "cubic 
spline  ".  (See  Ahlberg  et  al  section  2. 1 ).  We  have  to  choose  some  parameters  which 
characterize  the  polynomials  in  each  interval.  Let  these  parameters  be  s(x^)  =y^ 
and  s'(xj)=mj,  say.  This  is  a  meaningful  choice  since  s(x)  and  s’(xj)  are 
the  same  for  the  two  intervals  x^_j  £  x  £  x^  and  x^  s  x  £  xj  +  j.  Besides  that, 
the  four  quantities  yj,  nq,  yj  +  mj  +  ^  determine  the  four  parameters  of  a  cubic  polynomial 
in  £  x  s  +  1  completely.  The  continuity  requirement  for  *he  second  derivative 
at  points  x^,  x^,  •  •  •»xn_  ^  leads  to  the  following  set  of  equations. 

{  mi-i  +  2mi  +  imi+i  =  2ii(yi+i_yi-i)  <A-4> 


A  detailed  derivation  of  this  equation  is  given  in  the  above  stated  reference. 
However,  the  principle  of  how  these  equations  are  derived  should  be  clear.  The  two 
polynomials  in  x^_  j  £  x  <  x^  and  x^  £  x  sx^  +  1  are  linearly  expressed  in  terms 
of  y^f  Yj>  Yj+ji  mi_|»  mi»  mi  +  i  and  then  the  second  derivatives  at  x  =x^  are 
equated . 
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Since  there  are  in  (A.. 4)  n-2  independent  equations  for  n  unknowns,  two  more 
independent  conditions  are  needed.  They  are  chosen  involving  only  the  "boundary 
values"  nij,  m2  and  mn_  ptnn  .  However  we  shall  not  elaborate  on  this. 

Having  chosen  these  additional  conditions  in  some  way  the  quantities  mj  can 
be  solved  for  which  leads  to  expressions  like 


(A.5) 


It  is  now  important  that,  as  it  is  rigorously  proved  in  Ahlberg  et  al  (1967)  , 
section  2.4,  for  n  -*  <*>  and  i  -  n-i  -  ®,  (i.e.  for  a  large  number  of  points  and  i 
far  enough  from  the  boundary)  the  coefficients  tend  toward  a  limit  which 

depends  only  on  j-i  and  equals 


-  ai"1 


+  ai-j 


for  j-i>0 
for  j-i  =  0 
for  j-i  <  0 


(A.6) 


where 


=  -(2  -  V3)  -0.268 


(A. 7) 


Since  a  <  1  we  see  that  only  values  y\  at  in  the  neighborhood  of  x^ 
contribute  significantly  toward  m- ,  m^+  ^  and  therefore  to  the  polynomial  the 
interval  x.  £  x  s  x.,,. 

i  i  + 1 

(A.5),  (A.6),  (A. 7)  may  be  comprised  to 

\  ^<yi+j -yi-,>  <*•»> 

J  =1 

The  summation  limit  )3  is  chosen  in  a  way  that  is  sufficiently  small.  In  the  pract¬ 
ical  applications  of  cubic  splines  within  this  report  we  found  that  in  view  of  the  limited 
accuracy  of  the  data  (i.e.  heights,  gravity  anomalies)  a  value  of  j3  =  6  was  sufficient. 
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A. 2  The  Bi-Cubic  Spline 


For  a  function  depending  on  two  variables  a  generalization  is  necessary.  A  use¬ 
ful  approach  is  to  assume  that  function  values  are  prescribed  for  each  point 
(x^,  yj  )  of  a  rectangular  grid.  For  simplicity  we  assume  in  the  discussion  here  that 
the  spacings  x^  +  j-x^  and  y .  +  j-y^  are  equal  and  equal  to  h. 

One  can  proceed  in  the  following  way.  Cubic  splines  Sj(x  )  are  computed 
interpolating  the  given  z^  for  fixed  j.  Then  for  a  fixed  x  a  cubic  spline  s(x,y) 
viewed  as  a  function  of  y  is  computed  interpolating  the  Sj(x)  values,  in  that  way 
for  each  x,y  an  interpolated  value,  namely,  s(x,y)  is  obtained.  In  this  procedure 
the  role  of  x  and  y  (or  i  and  j)  may  be  interchanged  and,  surprisingly,  the  re - 
suit  is  the  same,  at  least  for  a  certain  class  of  boundary  conditions.  The  resulting 
function  s(x,y)  is  a  bi-cubic  polynomial  in  each  quadrangle  X|  :£  x  £  x^+p  £ 


y  s  yi+i»  i-e-: 


3  ,  -  -  v  .  x.-sx  sx.  , . 

s(x,y)  =  2  ak’Je>(x-x.)k(y-yj  )£  1 

k,  £  =  0  J  yj*y*yj+1 


(A.  9) 


s(x,y )  is  together  with  its  first  and  second  derivatives  continuous  over  the  entire 
range.  One  way  to  compute  the  is  the  following.  Perform  an  ordinary  spline 

interpolation  along  the  lines  y  =  y^ ,  x  =  x^ .  Obtain  from  that  according  to  section 
A. 2  the  values 


m. .  •- 

xj 


ds(x,yj) 


(A.  10) 


5s  (xpy) 
3y 


(A.  11) 


Perform  an  ordinary  spline  interpolation  of  m.j  for  fixed  i  and  obtain  from  that 


Pij  = 


5  s(x,y ) 
3  x  3  y 


(A.  12) 


x=xi ,  y  =yi 
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The  same  result  would  be  obtained  if  a  spline  interpolation  on  n^  for  fixed  j  was 
performed.  Now  compute  the  matrix  a1^  from  the  formulas 


a(i’J>  =  AT(h)K<‘-i>A(h) 


(A.  13) 


with 


A(h)  = 


1 

0 

0 


0  -3/m 

1  -2/h 


2/h3 

1/h2 


0 


3/h2  -2/h3 


0  0  -1/h  l/bf 


(A.  14) 


and 


,(i»j) 


Zi.j 

ni.j 

Zi.j+1 

Di»  j+1 

mi.j 

Pi,j 

mi»  j+1 

Pi,j+1 

zi+l,j 

Vi,j 

zi+l,j+l 

Di+l,j+l 

m,  ,  .  . 
l+l  ,j 

Pi+l.j 

mi+l,j+l 

Pi+l,j+l 

(A.  15) 


These  formulas  are  found  in  DeBoor  (1962)  and  have  also  been  used  in  Davis 
and  Kontis  (1970). 

It  is  clear  that  boundary  conditions  influence  the  values  m^,  n^,  p^  and 
therefore  a^  .  However,  this  influence  is  nearly  zero  if  from  now  on  we  assume 
that  the  grid  is  very  large  and  the  quadrangle  is  at  a  sufficient  distance  from  the 
boundary.  In  that  case  the  quantitives  m^,  n^,  p^  can  be  in  agreement  with  sec¬ 
tion  (A.  1)  and  can  be  computed  from 

m„  «  -  £  E 


nij  =  ' 


ij 


^  W 

(A.  16) 

(A.  17) 

(A.  18) 
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f 


or  alternatively 


P.. 

ij 


3 

h 


P 

£ 

k  =  1 


n.  .  ■) 

i-k.J7 


(A.  19) 


This,  together  with  (A.  13)  to  (A.  15)  was  the  set  of  formulas  used  in  this 
report.  Generalizations  to  grids  with  unequal  spacings  are  possible.  However,  one 
can  use  the  equal  spacing  formulas  also  for  an  originally  unequally  spaced  grid,  or 
even  for  a  curvi-linear  grid  applying  to  some  curved  surface.  The  derived  interpola¬ 
tion  polynomial  (A. 9)  has  only  to  be  interpreted  properly.  We  illustrate  tnis  in  the 
case  of  the  unit  sphere  with  a  0,  X  coordinate  system.  Let  the  grid  be  formed  of 
0,  X  lines  at  distances  of,  let’s  say,  1°.  This  grid  is  not  plane  nor  is  it  equally 
spaced.  Let  nevertheless  refer  to  the  grid  points  and  derive  the  polynomial 
(A. 9)  by  using  the  set  of  formulas  (A.  10)  -  (A.  15)  with  h  equal  to  tt/180  (arc 
length  of  one  degree  at  the  equator).  Let  s(9,X)  be  this  polynomial  where  €^X  (arc 
lengths)  replace  formally  x,y: 

3  3  (ij)  k  i  0ise*9i+l 

s(9,X)  =  E  £  a[lf,(d-Q)  (X-Xj)  (A. 20) 

k  =0  f  =0  **  Xj  ‘  X  ‘  Xj+1 


There  is  nothing  against  evaluating  s(9,X)  for  any  4X  value  not  coinciding  with  a  grid 
point  from  this  formula.  Also  derivatives  with  respect  t»  0,X  can  be  evaluated  by 
formally  differentiating  (A. 20).  These  derivatives  can  be  used  in  various  vector 
analytical  expressions  like  e.g. 


Lap  s(9,X) 


£s(£X)  + 

0  9  2 


MML  cot  0  + 

09  0X2  Sln20 


.  .  .  (A. 21) 
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APPENDIX  B:  Isotropic  Integral-Operators  With  Truncated  Kernel  Functions 


B. 1  Introduction 

Let  T  be  the  unit  sphere  and  let  the  unit  vectors  77  be  points  on  T.  A  function 
K(  5,77 )  is  called  a  kernel  function  since  it  may  occur  in  an  integral  transformation 


g(5)  =  jKts.ijMtoJdroi)  (B.l) 

r 


If  K(§,rj )  depends  only  on  the  distance  between  §,77  or,  equivalently,  only 
on  the  inner  product 


cos  0  =  5 . 17 


(B.2) 


then  K(£,  77)  will  be  written  as  K(  g-ry)  or  K(cos  0)  and  will  be  called  an  isotropic 
kernel.  It  is  known  that  the  spherical  harmonics  Snm(5)  are  the  eigen  functions  of 
such  integral  transformations: 


with 


n 


+  1 

2 IT  J  K(t)Pn(t)dt 
-  1 


<B.3) 

(B.4) 


PQ(t)  is  the  Legendre  polynomial.  (B.  3)  together  with  (B.  4)  is  called 

Funk-Hecke  formula  in  Muller  (1966).  If  the  spherical  harmonics  expansions  of 
f(£)  and  g(£)  are 


«!>  “  „Em  WW®  • 

n>m  n,m 


Then  (B.  1)  can  be  replaced  by 


g 


nm 


^n  nm 


(B.4a) 


Tliis  follows  immediately  from  (B.  1),  (B.  3).  Truncation  of  an  integral  transforma¬ 

tion  (B.3)  means  that  in  the  formula 

g(S)  =  /  K(§*  77)  f(r7 )  dT (77)  (B.5) 

r 
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the  integration  is  carried  out  only  over  some  (usually  circular)  cap  around  Let 
C  =  C(|,  0q)  denote  a  circular  cap  centered  at  |  and  having  spherical  radius  $q. 
Then  instead  of  (B.  5)  *.ve  have 

8 (?)  =  X  K(5*  t?)  f(r?)dr  (77)  (B.6) 

c 

Introducing  the  truncated  kernel 

f  K(5-tj)  for  §  *tj  a  cos  * 

K(?-tj)  =  0  <B.7) 

t  0  for  5*17  <  cos 

the  equation  (B.7)  can  also  be  written  as 

g(9  =  JK  (5*7?)  f(r?)  dT  (77) 

r 

This  shows  that  K (|  rj)  is  also  an  isotropic  operator,  hence  according  to  the  Funk- 
Ilecke  formula 

J  Z(Z-n)snm(r,)dr(V)  =  xnsnm(5) 

r 

a  _ 

xn  =  2v  r  K(t)Pn(t)dt 

-i 

The  error  function 

Ag(§)  =  g(0  -  8(5)  (B-10) 

is  then  given  by 

Ag(?)  =  J*  K( 5*  r?)  f (7? )  d T (77 )  =  J*  AK(5.r?)f(r))dr(rj)  (B.  11) 

r-c  r 

with  ~ 

A K(5*7j )  =  K(5-r))  -  K(§.rj)  (B*  12) 

Its  eigen  values  are 

AXn  =  Xn  ■  xn 


(B.8) 

(B.9) 
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with 


+  1 

AXn  =  2tr  J*  AK(t)Pn(t)dt 
-  1 


cos^ 

2v  /  K(t)Pn(t)dt 
-  1 


(B.  13) 


In  case  of  Stokes'  operator,  i.e. 

K(5-n>  =  -Lstu-n ) 

4tr 


with  Stokes  *  function  defined  by 


(B.14) 


the  AXr  are  given  by  ^  Qn  where  Qn  are  the  Molodensky  coefficients 

COS  ij) q 

Qn  ‘  i  St(t)Pn(t)dt 
-  1 


(B.  15) 


The  usual  method  to  evaluate  these  coefficients  is  based  on  Molodensky  et  a 1 
(1962)  and  is  designed  to  yield  Qn  =Qd($q)  for  fixed  n  and  all  ib^  .  In  section  B3-3 
we  shall  give  a  recursion  formula  which  will  yield  Qn($o)  for  fixed  j/jq  and  theoret¬ 
ically  all  n. 

We  shall  also  deal  with  the  following  question:  Is  there  a  way  to  decrease  the 
truncation  error  by  truncating  not  abruptly  at  ipg  but  continuously  over  an  integral 
ib  y  ^  ib  5  •  In  that  case  we  would  put 

K(cos$)  =  K(cos  0)  0  s  0  =s 

K(cos  0  =  0  0q  s  0  s  tr 

and  we  would  choose  K  (cos  0  somehow  in  the  interval  £0  ^  so  that 
K(cos  ib)  or,  equivalently,  AK(cos  0  has  some  desirable  properties  like  continuity 
or  continuous  differentiability. 

It  is  known  that  generally  the  Legendre  coefficients  decrease  more  rapidly 
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to  zero  as  n  goes  to  infinity  if  the  relevant  function  is  smooth.  Therefore  AAn,  the 
eigen  values  of  the  error  kernel 

AK(cos$)  =  K(cos0)  -  K(cos  $)  (B.16) 

can  be  expected  to  taper  off  more  quickly.  This  is  for  example  seen  from 
+ 1 

AAn  =  2f r  J  AK(t)Pn(t  )  dt  =  J  AK(5.n) 

‘  1  r 

by  Green’s  theorem  (2.5)  for  n  si: 

=  J  Lap  AK(  5  -r))  Lap  P  (?*7j)dr  (n)  = 

0  rj  n  11 

r 

-  1  J*  Lap  AK(5-r})  Pn(|.rj)  d  r  (rj) 

n(n  +  l)  p 

If  AK  is  two  times  differentiable  so  that  the  surface -Laplacean  Lap  AK 
exists  and  is  squared  integrable,  then  it  follows  from  the  Schwarz  inequality  that 

i - j - - - > 

AX"  s  i^Ti)  y  j „» 2dr<>))y;(pn(5^))  »r(,) 

Thus,  since  the  integral  involving  Pn  is  O(-i-) : 

A\n  =  0(__L_)  (B.  17) 

n2Vn 

In  section  (S.4)  we  shall  deal  with  this  question  more  systematically.  One  of 
the  necessary  preparations  will  be  to  compute  the  surface  Laplacean  for  certain  fre¬ 
quently  used  kernel  functions  K(5-  77). 


13 . 2  Tie  Surface  Laplacean  of  Some  Kernel  Functions 

If  K(^rj)  is  an  isotropic  kernel  function,  we  are  interested  in  the  quantity 
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K(5*i?) 

which  denotes  the  surface  Laplace#  n  with  respect  to  77.  (Because  of  the  complete  sym- 
metry  of  the  problem  the  Laplacean  with  respect  to  5  would  give  the  same  result.) 

If  we  put 

T 

17  =  (sin  9  cos  X ,  sin  9  sin  X ,  cos  X ) 

and  if  F(9,X)  denotes  a  function  on  r,  depending  now  on  polar  distance  angle  0 
and  longitude  X,  then  the  Laplacean  may  be  computed  from 

Lap  F<e.  X)  =  F00  +  F0  cot  9  +  <»•»> 

In  our  application  we  assume  the  north  pole  of  die  6^X  system  in  %.  K(|*  rj )  may 
then  be  written  as  K(cos  0)  and  is  independent  of  X.  Thus 


=  Lap  K(cos  0)  =  +  KQ  cot  0  (B.  19) 

B.2-1  The  Kernel  l/HZ-n).  Let  4(?«fj)  denote  the  distance  between  %.r), 
i.e.  with  5  -rj  =  cos  <p 

4(5*1?)  =  5-1?  =  2  Sin  £  (B.20) 


Straight  forward  differentiation  according  to  (B.  19)  and  (B.20)  with  0  =  0  yields 
the  result 

1 _  =  1  4-  1 


Lap 


V  4(§.t?)  f3(§.r?) 


44(5*1?) 


(B.21) 


B.2-2  The  Kernel  4(5*  rr). 

Lap 

i? 


In  an  analogous  manner  we  obtain 


4(5*1?)  = 


4(5*r?) 


MLXll 


(B.22) 
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9  *  2n+l 

4  n  =2  (n  +  i.)n(n  +  l) 


Pn(?*r?) 


+  3  Z 


2n  +  1 


n  =  2  (n  -  1 )  (n +^)n(n+l) 


pn<5-n) 


_  (B.28) 
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which  may  be  verified  taking  into  account  the  spherical  harmonics  expansions  of  the 
involved  functions.  The  idea  behind  (B.28)  is  that  the  La  place  an  can  be  applied  term* 
wise  on  the  right-hand  side  while  the  left-hand  side  causes  no  problems  in  view  of 
(B.21)  and  (B.23).  Performing  the  necessary  elementary  manipulations  one  arrives 
finally  at 

l^P  St(S-r?>=  -5^ - 2St(?*7>)  +  2P  (§-T7>  +  9P1(?*t?)  (B.29) 

£3<?*77)  ° 

Remark:  The  procedure  can  easily  be  extended  to  residual  Stokes'  kernel  functions 
which  are  obtained  after  eliminating  from  St(^*tj)  spherical  harmonics  up  to  and  in¬ 
cluding  degree  N.  Call 

StN(5-7j)  =  StU-rj)  -  E  ^±-1Pn(?-rj)=  E  2J?.±A,  P  (g-rp)  (B.30) 

n -2  11-1  n=N+ 1  n_1  n 


Then  0?.  29)  generalizes  to 


N 


Lap  St  (§-rj)  =  - 2StNr(?.T]i>  +  E  (2n2  +  5n  +  2 )  P  (l-r)  (B.31) 

N  r(5*t|)  n=0 


B.3  Truncation  at  a  Certain  Angle  i/j^. 

In  this  section  we  concentrate  on  the  error  kernel  AK(|*  r\)  defined  in  (B.  16) 
having  eigenvalues  AXn  according  to  (B.  13): 

+ 1  cos  jL 

AX n  =  2 tr  J*  AK(t)Pn(t)dt  =  2ir  J*  K(t)Pn(t)dt  (B.32) 

-  1  -  1 

If  the  truncated  kernel  is  applied  to  a  function  f(§  )  ~f  (f^  -  spherical 

harmonics  coefficients)  then  the  resulting  error  Ag(5)~Ag  is  given  by 

nm 

Ag(?)  =  J*  AK(?*tj)  f<77 )  d  T  (rp)  (B.33) 

r 
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or  in  view  of  (B.4a)  by 


Ag  =  AX  f 
^nm  n  nm 


(B.34) 


We  shall  now  compute  the  AX  's  for  some  kernels: 


n 


B.3-1.  The  Kernel  1/j g(g-n):  From  (B.20)  and  (B.32)  we  get  with  cos  ($q)  =  t^ 
because  of  sin  =\J  -  ^os  ^ 


AXn  -V5»  J  —4=  Pt(t)dt 
- 1  V  i  - 1 


(B.35) 


Write 


^n  X 


m\  vt--i  n 


—  P„  (t)dt 


(B.36) 


Then 


»  bn 


(B.36a) 


We  derive  now  a  recursion  formula  for  bn-  First,  for  n  =  0,  1  one  computes 
in  an  elementary  way 


b„  *  -  2  Vl  '  «(,  +  2VT 


b 

1 


(B-37) 


For  n  2  2  we  replace  Pn( t )  by 


p„">  - 


n-1 


Pn-2W 


cf.  Abramowitz-Stegun  (1964),  formula  (8.5.3).  Inserting  into  (B.36)  yields 
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b  = 
n 


2n-l 

n 


c0 

J 


1  Pn_l(t)dt  -  — -  b 


v~r^ 


n  n  -  2 


Ir  S  VTr-t  pn.l(t)  „t  +  2fib 

- 1  n 


t 0 


Applying  partial  integration  and  using 
♦1 


I  P„<'> 

-1 


dt  = 


2n  +  l  '  n  +  rlO'  n-1'  0 


(cf.  Lense  (1954),  p.  17)  we  get 


b  _  Pn-2  «0>  ~  Pn  <«0> 
"  n 


•  27<b„  •  bn-2>  +  ^Vl 


solving  for  bfl  yields 


n 


Pn-2(t)  ‘  Pn(t)  W -  **  2 

b-  =  — - - -  V  1  "  c0  +  2 - ^~b 


n-1 


n  +  JL 
2 


n  bn-2 


3 

n‘Tb. 


n+  1  n-1  n  ,  1  ”n-2 

n+  -  "  +  j 


(B.38) 


This  is  the  desired  recurrence  relation.  bn  could  be  replaced  by  using  (B.  36a) 

B.3-2  The  Kernel  l/l^g.^).  For  this  kernel  we  have 

tr 


L0 


AXn  ~  \p>  «f 


1 


V2  «  \/  ( 1  -  t)d  n 

-  1 


3  PJO  dt 


(B.39) 


Writing 


AA_  = 


IT 


n  "  VT  n 


(B.40) 


with 


*  *  VTT7?)3  P"<t)dt 


(B.41) 
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One  derives  in  an  analogous 


way  the  recurrence  relations 


C0  = 


(B.42) 


:i  = 


i-  2  Vi  -  tr 


=  Pn-^Q)  -  Pn(tn) 


+  2  c 


Stokes'  Kernel,  y/e  take 
(B.32)  in  the  form 


a  slightly  different  approach  and  write  the 


0.43) 


equation 


AXn  -  J  st(§. 77)  Pn(S-T])<ir  (rj) 


r-  C  is  d,e  complementary  area  after  removal  a  cap  with  spherical  radius  #  from 
the  unit  sphere  r.  Nert  we  apply  Green’s  second  formula  and  obtain  for  n  n  1  ° 


Mn  ~  1  ^p^lpn(;-n)drip) 


§St(%Tj,  (Omd^Lap'1  „<„>  >d/S(„> 

P 

§  Ulfv  Pn(5->l><Grad^  St<;-„),  ulnDtiu  (r?) 


f  1S  the  boundary  of  r-C  and  p  is  the  unit  normal  ,0  the  boundary,  tangen¬ 
tial  t„  rand  outward  of  r-C  (inward  c).  Choostng  a  0,  X  coordinate  system  with 
P0.e.n  (,  snbsuntung  m  the  line  Integrals  t  -  cos  9,  noth*  -(1**  +  1)) 

n  “  *1'  ’  intr°ducing  tQs  cos  cp^t  and  verifying 

(Grad_  F(  jj),  v(r>))  =  -  =  4.  d  F(t)  \  H  2 

77  d  e  iu  V 1  ‘  f 
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leads  to 

AXn  =  •  ToT+T)  f  Lap^std-^p^^TjJdr  (n) 

r-c 


2  ff 

n(n  +  l) 


St(tQ) 


dPn(t) 

dt 


(1 


t=tr 


l0> 


2 ff  d  St(g-  n) 
n(n  +  l)  dt 


t=L 


PJtrt)(l-tft) 


n 


(B.44a) 


Denoting  the  derivatives  with  respect  to  t  by  dashes  using  (B.29)  for  Lap  St 


gives 

AX 

n 


T*rn>  I  ^3—  pn<5--7)«r(,)  +52—  X 

r-c  r  -  c 

^f+T)  f  p„(?-’7>dI'<n>  -  f  P,(s-»>P.(5-n>drW 

r-c  r"c 

^l+T)  St<'o>  Pn *t0*  < 1  ■  'o*  +  ^T)** '<•«>>  (B-44> 


Using  the  results  of  section  B.3-2  and  equation  (B.44)  leads  to 


AXn.=  " 


£rr» c-  +  AX” '  »o>+5  p-<t)dt 


t  =  l 


18  7T 

n(n  +  l) 


f  tPn(t)dt  -  -42rr.St(^)P,(t„)(l-^) 


t=“  1 


-^71)  vlo' vLo' 


^T)*t,(i0>Pn<*0>(‘-«S> 


Replacing  AX  by  the  Moiodensky  coefficients  Q  according  to 
n  n 

AXn  =  2ffQn  (B.45) 
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and  solving  for  Qn  gives  finally 


Q 


n 


1 

n(n  + 1 )  -  2 


c 

n 


n(n  + 1 )  -  2 


1  2  ^ 

Z  <2mz+5m+2)  J*  Pm(t)  Pn(t)  dt 

m  =0  -  1 


(B.46) 


n( 


St<'o>  p:<'0>  <*  -  *0>  +  -„7-n%7)-  2  st  '<‘0>  V‘o>  <*•«;> 


This  is  the  desired  formula.  (It  holds  also  for  StN(?*rj)  instead  of  St(§-r^ 
if  the  summation  over  m  is  extended  to  N.)  See  (B.  30)  ,  (B.31).  Formula  (B.46) 
is  recursive  insofar  as  the  cn  are  the  result  of  the  recursion  (B.42)  -  (B.43).  The 
two  integrals  involving  Pn(t)  may  be  computed  from  the  general  formula 


- 1 


Pm('0)  t 


n«0>  -  WPn<*0> 


n(n+l)  — 


m(m  + 1 ) 


for  n/m 


(B.47) 


This  formula  follows  from 

« 

+  1 

f  V'>  V‘>dt  -  it"  J-  pm<5-”)pn<?^dr(n) 

- 1  r 

after  application  of  Green’s  second  formula.  It  holds  for  m  ^  n,  i.e.  it  may  be  used 
for  n  >  2  in  (B.46)  for  n  =  0, 1  elementary  integrations  are  performed.  The  deriva¬ 
tive  PjJ(t())  is  obtained  from 

p'to  =  -  -^(tPn(t)  -  Pn  .(t))  0.48) 

n  j  _  j.  2  n  n- 1 


(Lense  (1954),  p.  16)  or  a  similar  formula.  The  derivative  of  Stokes'  function  is  in  a 
direct  way  computed  from  (B.24)  observing 


St'(t)  =  d  St(9°8.’P) 

dt!) 


-1 

sin  il) 


cos  tj) 


t 
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0.49) 


B.4  Modified  Truncation  Procedures.  Repeating  formula  (B. 44a)  for  a  general 
kernel  K(§*77)  we  get: 


* 


'  ^i+T)(l’tO)KitO)P'nitO)  +  ^T)<1“t0>K,<t0>Pn<t0>  <B*50> 

One  sees  that  the  first  and  also  the  third  term  on  the  right-hand  side  taper  off  like  an 
0(— Jj-)  as  n-®,  whereas  the  second  term  which  contains  P^(tQ)  does  not  necessarily 
in  view  of  (B.48).  It  is,  therefore,  natural  to  try  to  modify  the  kernel  K(cos  (p)  in 
a  way  that  K(tg)  vanishes.  This  can  be  accomplished  by  simply  subtracting 


kQ  =  ICO*,) 


The  effect  on  the  original  integral 


g<5>  =  jK(§^)f(rJ)dr(r?) 


(B.51) 


(B.52) 


is  independent  of  5  and  therefore  constant.  For 


g(s>  =  J*(k (5-t))  -  k0)f(n)dr(n)  =g(§)-r0  (b.53) 

r 


70  =  ko  /  f(rj )  d  T  (77)  =  4trkQf 


(B.54) 


f  denotes  the  mean  value  of  f (17 )  over  T. 

In  most  applications  a  constant  error  in  g(£)  does  not  matter  at  all.  to  Stokes’ 
formula  for  example  the  harmonic  of  degree  zero  is  eliminated.  The  scale  of  the  com¬ 
puted  geoid  has  to  be  obtained  from  other  sources.  Therefore  a  constant  error  is  with¬ 
out  concern.  It  is,  however,  important  that  the  truncation  angle  is  the  same 
everywhere . 
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If  K(t)  happens  to  have  a  zero  at  t  =  tg,  then  no  subtraction  is  necessary. 

De Witte  (1966)  found  small  truncation  errors  at  the  zeros  of  Stokes'  function.  The 
theoretical  reason  for  this  is  clear  by  now.  It  is,  however,  also  clear  that  the  error 
at  any  other  truncation  angle  is  to  a  large  extent  a  constant  over  all  of  r  if  kg  is 
subtracted  prior  to  truncation. 

One  can,  of  course,  also  try  to  modify  die  kernel  K(cos  if)  in  a  way  that  also 
the  derivative  at  the  angle  of  truncation  vanishes .  Let  this  angle  be  again  denoted  by 
j|)q.  Then  K(cos  if))  will  be  modified  in  s  if)  s  $  so  that 

K(cos  if))  =  K(cos  if))  ,  if)  s  0^ 

K(cos  t/j)  =  0  ,  if)  z  if) q 

In  ip^  <:  if)  ^  if) q  K(cos  ip)  will  be  chosen  in  a  way  that  K(cos  if>)  is  twice 
differentiable  in  0  <  <p  z  tt.  The  error  kernel  is  then 

AK(cos  if))  -  0  ,  if)  £  il)^ 

AK(cos  if)  =  K(cos  if)),  if)  2 

and  AK(cos  <p)  twice  differentiable  in  0  s  ^  s  n. 

The  question  is  now,  how  to  choose  K  or  AK  within  the  interval.  Since  the 
eigen  values  of  the  error  kernel  are  now  by  (3. SO)  : 

=  '  -^ij  J^^(5^)Pn(5  n)dr(,) 

r  (B.55) 

+  1 

=  '  nT^Tn  S  Lap„AK(t)P  (t  )dt, 

n(n  +  i)  «  r)  n  J 

-  1 

it  will  be  beneficial  if  the  Laplacian  of  AK  is  as  small  as  )x>ssible.  Since  the 
Laplacian  depends  on  the  second  and  first  derivative  it  may  be  good  strategy  to  try 
to  make  these  derivatives  small. 
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We  view  the  kernel  AK  as  a  function  of  t  =  cos  if).  With  tQ  =  cos  <pQ,  tj  =  cos  <pv 
we  take  a  piecewise  quadratic  polynomial  AK(t )  in  tQ  <  t  s  tj  such  that 

AK(t())  =  K(t0)  ,  AK  '(Iq)  =  K'(t0) 


AK(tp  =  0  ,  AK*(tj)  =  0 


(B.56) 


and  AK  (t )  is  a  polynomial  of  degree  2  in  each  of  the  intervals  tQ  st  <;  t*, 
i*  t  i*  will  be  chosen  so  that  AK11  is  of  equal  modulus  (but  opposite 

sign)  in  the  two  intervals.  Of  course  AK(t  )  and  AKf(t)  shall  be  continuous  at 

t  =  t*. 


We  do  not  write  down  the  explicit  set  of  formulas  which  is  completely  elementary. 
It  may  only  be  mentioned  that  t*  may  coincide  with  t^  which  appears  desirable  since 
then  one  of  the  discontinuities  of  the  second  derivative  of  AK(t)  disappears.  This 
special  case  may  for  given  tQ,  t^  be  obtained  by  adding  a  constant  to  AK(t).  Adding 
a  constant  to  AK(t  )  amounts  to  adding  one  to  K(t  )  which  is  of  not  much  concern  as 
we  have  seen  earlier. 

More  sophisticated  kernel  -  modifications  can  be  conceived.  For  example  the 
second  and  third  derivatives  of  the  truncated  kernel  (or,  equivalently,  of  the  error 
kernel)  could  also  be  made  continuous .  This,  however,  will  not  be  done  here. 
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where  ct(?,  tj)  denotes  a  unit  vector,  tangential  to  T  in  the  point  ?  and  being  in  the 
plane  spanned  by  ?  and  rj •  Formula  (B.  62)  is  the  result  of  the  differentiation  of  a 
parameter  integral  where  the  parameter  §  occurs  also  in  the  domain. 


Let  Xn  denote  the  eigen  values  of  the  truncated  kernel.  Cf.  (B.8-9).  Then 
(B .  6 1)  may  be  written 


g(i) 


=  E 
n,m 


g  S  (?) 
nm  nm 


£  X 
n,  m  n 


f  S 
nm  nm 


a) 


(B.64) 


From  (B.60)  follows  then 

w(?)  =  E  g  GradS  (?)  =  E  X  f  _  Grad  S  (?)  (B.65) 

nm  nm  _  _  n  ntn  nm 


This  means  that  g  =  X  f  are  the  coefficients  of  w(  ?)  with  respect  to  the  system 
6nm  n  nm 

GradS  (?).  This  system  is  orthogonal.  Cf.  Meissl  (1971),  section  2.  It  is,  how- 
nm 

ever,  not  normalized,  even  if  the  S  are.  If  the  S  are  normalized  in  the  usual 

nm  nm 

way,  then  the  mean  square  norm  of  g(?)  is  given  by 

||g||2  =  J*  g2(5)d  T  (?)  =  4ff  E  g2  =  4tt  E  (Xnfnm)2  (B.66) 

r  n,m  11111  n,m 

and  that  of  w(  ?)  is  given  by 

II  wl!  2  =  J  (w(?),  w(?))dr(?)  =  4*  E  n<n  +  1>^nfnm>2  0.67) 

„  »,m 


This  is  more  fully  explained  in  the  stated  reference. 

Because  of  the  linearity  structure  the  above  discusrion  carries  immediately 
over  to  the  truncation  errors . 

Ag(!)  =g(5)“g(l)  »  Aw(?)  =  w(?)-w(?) 

Let  AXn  denote  the  eigen-values  of  the  error  kernel  as  in  (B.  13).  Then 
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(B.64a) 


Ag(?)  =  £  Ag 


nm  nm 


(?)  = 


AX  f  S  (?) 
n  nm  nm'5 


Aw(?)  =  E  Ag  Grad  S  (?)  =  £  AX.  f  Grad  S  (?)  (B.65a) 

n,m  nm  11111  n, m  n  nm  nm 

||Ag||2=  J*  Ag2(?)dT(?)  =  4ir  E  Ag*  =4tr  E  (AXn  fflm)2  (B.66a) 

n,m  n,m 

||  Aw||  2  =  J  (Aw(? ),  Aw( ?))  d T  ( ? )  =  Alt  £  n(n  +  l)  (AXn  f  )2  (B.67a) 

"  n,  m 


The  relationship  (B.65)  is  quite  simple  and  offers  an  easy  dicussion  of  the  trun¬ 
cation  errors.  The  usual  truncation  of  the  Vening-Meinesz  integral  is,  however,  not 
based  on  (B.62),  but  rather  on 

1v(?)  *  j*  Grad^C(?.f?)f(t))dr(i))  (B.68) 

C? 

De  Witte  (1966)  has  noted  this  difference.  He  used  associated  Legendre  -functions 
to  deal  with  formula  (B.68).  hi  our  treatment  here  a  quite  simple  analytical  expression 
for  the  difference  t»e tween  (B.68)  and  (B.62)  is  at  hand: 


w(?)  =w(?)  -  J  K(?-tj)  a(?,T7)f(r?)  d  a  C^(r7) 
ac„ 


(B.69) 


This  shows  for  the  first  that  the  difference  vanishes  if  K(cos  0)  happens  to  have 
zero  at  0  =  0^  .  This  is  quite  apparent  from  the  graphs  in  De  Witte  (1966),  though  he 
was  seemingly  unaware  of  the  simple  relationship  between  the  two  different  truncation 
procedures . 

We  tc-n,  of  course,  enforce  a  zero  of  the  kernel  at  0q  by  simply  subtracting 
the  constant  kQ  =  K(cos  </^)  as  in  the  discussion  around  (B.51).  The  effect  of  this  is 
a  constant  in  (B.59)  and  zero  in  (B.60).  A  beneficial  effect  is  noted  if  the  formulas 
are  truncated.  As  we  have  shown,  the  AXni.e.  the  eigen  values  of  the  error  kernel 
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taper  off  more  quickly.  Moreover,  the  two  quantities  w(§)  and  w(§)  coincide.  This 
shows  that  the  truncation  w  can  be  considered  as  favorable  too,  at  least  for  the  higher 
harmonics.  Equation  (B.67a)  shows  that  the  higher  harmonic  contribution  to  Aw  = 

Aw  is  in  the  same  dampened  by  the  factor  AX  as  in  the  case  of  Ag. 

Li  we  deal  with  Stokes  and  Vening-Meinesz  formulas  then  some  constant  factors 
enter  the  discussion.  For  we  have 

Stokes:  N<?>  =  ^  J*  St<5  -rj) Ag(rj)  d  r  <r>)  0.70) 

r 

Vening-Meinesz : 

v(5)  =  "  J*  Grad-St(l-r))  Ag(rj)  d  F  (17)  0-71) 

r 

These  factors  are  easily  taken  care  of. 

Summarizing  we  may  say  the  following:  Adding  the  constant  -  St(cos  ^  )  to 
Stokes'  kernel  causes  a  constant  added  to  the  undulation  N(f )  which  is  of  no  concern 
since  the  geoid  has  to  be  scaled  later  on  anyway.  The  added  constant  drops  out  completely 
in  the  Vening-Meinesz  formula.  If  the  formulas  are  truncated  at  0^,  then  the  higher 
harmonics  of  the  truncation  error  are  favorably  dampened.  Moreover  the  usual  trunca¬ 
tion  of  the  Vening-Meinesz  formula  yields  exactly  the  slope  of  the  geoid  comjxited  from 
the  truncated  kernel  St(  £  •  rj)  -  St(  cos 

If  the  Stokes '  kernel  is  truncated  over  an  interval  s  0  £  Yq  so  that  the 
derivative  at  y  =  y  a*so  vanishes,  then  the  same  can  be  said.  The  dampening  effect 
may  be  a  little  more  favorable.  The  error  discussion  of  the  Vening-Meinesz  formula 
can  be  based  on  the  various  quantities  AXn  as  they  are  discussed  in  sections  B.4,  B.5. 

13. 6  Smoothing  the  Integrand  in  the  Outer  Zones.  Up  to  now  we  have  been  treating 
the  case  where  the  integrand  is  completely  neglected  outside  a  certain  area.  In  this 
section  we  discuss  a  refinement  in  the  sense  that  the  integrand  is  not  replaced  by  zero 
but  rather  by  a  smoothed  version  like  a  moving  average  over  a  certain  area. 


56 


Let  us  start  with  formula  (B.  1)  which  we  write  now  as 

g((V2)  =  J* K(0)( £ -7j)  f  ^0)(rj)  d  r (rj)  (B.72) 

r 

~(0) 

Assume  truncation  (in  some  way)  at  an  angle  0q.  Denote  by  K  the  truncated 

kernel  and  by  AK^(§.?])  the  error-or  residual  kernel.  Then: 

K(0)<5*t?)  =  K(0)(5-tj)  +  AK(0)(§-n>  (B.73) 

The  error  is  then 

Ag<0)(?)  =  J*  AK(0>(?.7,)  f<°>(r,)dr(rj) 

r 

or  according  to  section  B.  1: 

-  AX<0>f<°> 
nm  n  nm 

with 

AX^0)  =  J*  AK<0)(5*rj)  Pn  (I  n)  d  T(t7)  (B.76) 

r 

The  shape  of  the  error  kernel  depends  on  the  philosophy  of  the  truncation 
procedure.  Confer  sections  B.3  to  B.4. 

Let  us  account  for  the  truncation  error  by  computing  the  correction  term  (B.74). 
If  we  do  this  by  using  f^(f)  in  its  original  form  then  of  course  no  truncation  error 
remains.  We  could  have  used  the  non-truncated  kernel  in  the  first  place.  The  point 
is,  however,  that  in  the  correction  formula  (B.74)  a  simplified  version  of  f'^(5)  may 
be  used.  Let  us  use  f^(5)  instead  of  f^(§)  in  (B. 74)  and  call 

Af<°><?)  =  f(0)(§)-  f(i)(5)  (B.77) 

then  a  residual  error  of 

Ag(1>(?)  =  j*  AK(0)(§.r?)Af(0)(T?)dr(T7)  (B.78) 

r 


(B.74) 


(B.75) 
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remains.  The  spherical  harmonics  of  this  error  can  be  computed  by 


Ag<‘>  =  Axf  Af<°> 
nm  nm 


(B.79) 


Assume  that  f^(  ?)  is  the  result  of  the  application  of  an  isotropic  smoothing 
operator  toward  f^(§).  If  p^  are  the  eigenvalues  of  this  smoothing  operator, 
then  we  have 


fU)  _  fi(0)f(0) 

nm  n  nm 


(B.80) 


and  consequently 

Af<°>  .  <  !-/»>,«»  =  V0) 

nm  n  nm  n 

so  that  (B.79)  becomes 

Ag<‘>  =  AX®  //>!« 
nm  n  n  nm 

The  A  X<°>  generally  taper  off  with  increasing  n 
while  doing  this).  From  a  certain  n  on  they  will  be  negligible  under  any  circum¬ 
stances.  It  is  therefore  desirable  to  have  small  A/J^  for  moderately  large  n.  This 
means  one  should  have  p^  close  to  1  for  moderately  large  n.  One  way  is  to  put 

P^  =  1  for  certain  nsnn  and  p^  =  0  for  n>nn.  This  amounts  to  a  replacement 

n  (0)  n  (1) 
of  fv  (§)  by  its  spherical  harmonic  expansion  f'  '(5)  truncated  at  n=np. 

More  common  is  another  smoothed  version;  namely,  averages  over  certain  block 
areas.  A  theoretical  discussion  of  this  is  somewhat  complicated  since  this  smoothing 
operator  is;  1)  not  isotropic  and  2)  applied  only  at  discrete  locations  (which  may  be 
taken  at  the  mid  points  of  the  biocks).  Easier  is  a  discussion  based  on  moving  averages 
over  circular  caps  of  half  opening  angle  Fbr  error  consideration  it  is  justifiable 

to  replace  the  block  averages  by  moving  averages  over  circular  caps  of  comparable 
size .. 

The  eigen  values  /3^  of  these  operators  are  found  in  Meissl  (1971),  equ.  (3. 14a), 
(3.  14b).  See  also  Pellinen  (19(56).  We  use  the  above  cited  formula  (3. 14a): 


f(0) 


nm 


(B.81) 


(B.82) 


(though  they  may  oscillate 
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(B.83) 


?(0) 

n 


1 - f  P_(t)  d  1 

1  -  cos  a0  d  n' 


cos  a, 


0 


For  moderate  large  and  n  we  may  replace  PQ(t)  by  its  Taylor-linearization  at 
t  =  1: 


P  (t)  =  1  -  n<n.+  1)  (1-t)  +  0( ( 1-t )2) 

2 

This  leads  in  a  straight-forward  way  to 

0(0)  =  i  -  n(n  +  l)  (1  .  cos  +  0((1  -  cos  c^)2) 


(B.84) 


(B.85) 


Neglecting  the  0-term  we  have 

*  1  -  ( 1  -  cos  a0)  (B.86) 

For  the  purpose  of  error  estimates  this  formula  may  be  used  as  long  as 

/J(0)  >  0.80,  say. 
n 

Formula  (B.72)  is  now  replaced  by 
g(0)(?)  =  = 

r 

=  jK(0)(5-r))f(0)(n)d  r  (r7)  + 

r 

+  jAK(U)(^r7)f(1)(r7)d  r(rj)  +  Ag(1)(|)  (B.87) 

r 

with 

Ag(1)(?)  =  jAK(0>(5.T?)Af<%)d  r(r?)  (B.88) 

r 

'The  following  nota'  ons  have  been  adopted 
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original  kernel 


k Vhz-v) 

AK(°>(|.t7) 

f(°>(5) 

f(1)(?) 

Af(0)(P) 


truncated  kernel 
residual  kernel: 
original  function 
smoothed  version  of  function 
residual  f^(§)  -  f^(§) 


It  is  considered  that  the  first  two  terms  on  the  right-hand  side  in  (B.87)  are 
computed.  The  third  term  Ag^7(§>  represents  the  error  term  (B.88)  which  in  spher¬ 
ical  harmonics  representation  is  given  by  (B.82). 

One  can,  however,  go  a  step  further  and  try  to  truncate  AK^(§*tj)  again  further 

outside 

AK(0)(S.t7)  =  AK(0)(?.r?)  +  AK  (B.89) 

The  contribution  of  AK^(5’Tj)  could  be  taken  into  account  but  only  after  replacing 
f^(5)  by  a  still  heavier  smoothed  version  f^\§). 

There  is  no  need  for  a  detailed  discussion  of  the  further  procedure  which  could  go 
on  and  on,  replacing  f^(§)  by  smoother  and  smoother  versions  in  a  succession  of 
concentric  zones.  We  shall  restrict  us  here  the  following  remarks: 

(1) :  If  truncation  over  an  interval  is  used,  then  there  is  an  overlapping  of 

neighboring  zones.  This  does  not  matter  if  it  is  properly  taken  care  of. 

(2) :  One  should  avoid  adding  constants  toward  the  AK^  kernels  before  trun¬ 
cation  as  it  has  in  some  cases  been  done  for  the  original  kernel  K^(§.  77).  This  would 
make  the  kernels  non-zero  also  in  zones  farther  inside.  If  AK^  happens  to  have  a 
zero  at  a  location  which  is  suitable  for  truncation  then  truncation  is  advantageous  in 
the  same  way  as  it  was  for  a  zero  of  K(cos  0).  Kernels  from  which  the  harmonic 
components  up  to  degree  N  have  been  removed  have  at  least  N+l  zeros.  (Meissl 
(1971)  ).  These  zeros  are  natural  locations  for  truncation,  or  for  a  transition  to  a 
heavier  smoothed  version  fO  +  l)(F)  of  f^^(P),  respectively. 
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APPENDIX  C:  Regularization  of  a  Type  M- Integral 


What  we  call  a  type  M-integral  looks  like 

g<§>  =  j*  dr  (77)  (C.l) 

r  *  (Zri) 

It  is  singular  and  we  require  that 

<p(z,  0  =  0 

and  also  that  <p(§,  r?)  is  twice  continuously  differentiable  throughout. 

Examples  are 

g(5)  =  J  f^j  d  r  (77) 

0  jg'3 7 ^ 

r 

or 

g(§)  =  J  ^  h(T?)d  r (77) 

r  ^ 

The  latter  is  the  correction  term  in  the  first  order  solution  of  Molodensky's  problem. 
Therefore  the  name  "type-M"  has  been  chosen. 

The  regularization  procedure  will  be  based  on  Green's  second  formula.  Only 
a  cap  C  will  be  considered  since  for  larger  distances  the  integral  tapers  off  quickly. 
There  the  original  form  may  be  used.  Complete  truncation  at  a  certain  distance  is 
also  feasible  in  most  applications. 

The  contribution  of  a  spherical  cap  C  centered  at  £  and  having  half-opening 
angle  ^  to  a  type-M  integral  is 

gc(S>  =  S  d  T  (77)  (C.5) 

C  c£3U-r?> 

Eonmila  (B.21)  allows  us  to  rewrite  this  as 


(C.2) 


(C.3) 


(C.4) 
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(C.6) 


gc(5)-gj(5)+g2(S)  =  J*  Lap >v  —  §•  77)  d  T (77) 


*1 


1  1 


4  j e  ( §  *77) 


<p(§»r?)  d  r(7j) 


The  second  integral,  i.e.  that  one  for  g  (5)  is  already  regular.  In  some  ap- 

2  .  1 

plications  it  may  even  be  omitted  since  after  multiplication  by  R  it  is  negligible 
in  the  so-called  planar  approximation.  See  Moritz  (1969). 

Before  we  apply  Green's  second  formula  toward  g  (|)  we  exclude  a  small  cap 
C|  around  5  of  half  opening  angle  ij)^.  We  have  then 

gjCS)  =  J  Up»7(7^j  * 

c  -c,  ' 


■  s 

c-c 


■■  Lap-  vi)  d  r  (77)  + 

Ul-T})  * 


+  §  <p(Z>rj)  (Grad^  »  v(lj)  )  d  3  C(tj)  - 

ac 

-  |  — i —  (Grad  <p(  §,77),  v  (77)  )  d  9  C  (77 )  - 
*<S-T7>  V 


ac 


|  <P  (?» 17)  (Grad^  1^(77)  )  d  3  C  ^ij)  + 


a  c  1 


+  §  T, - 7  (Grad  77),  i/,(tj  ) )  d  a  C  (77 ) 

o  t- , 


(C.7) 


v(n  )  is  (outward)  normal  to  the  boundary  of  C  and  so  is  1^(77)  to  the  boundary 
of  r.  If  we  let  C  contract  toward  the  point  %  then  it  may  be  shown  that  the  last 

*  I 

two  integrals  vanish  because  of  the  assumptions  regarding  <p( 77) .  The  third  integral 
on  the  right-hand  side  may  be  combined  with  the  first  one.  This  is  possihle  since 
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APPENDIX  D:  Regularization  of  the  Vening-Meinesz  Formula 

We  transform  the  Vening-Meinesz  formula  in  a  way  that  the  resulting 
integrals  are  no  longer  singular.  We  deal  thereby  only  with  the  contribution  of  a  cir¬ 
cular  cap  C  around  the  point  5  in  which  the  deflection  is  to  be  computed.  The  hali- 
opening  angle  of  this  cap  shall  be  j/>^.  For  the  more  remote  zones  die  original  form 
of  Vening-Meinesz  formula  is  of  greater  advantage  since  the  kernel  tapers  off  quickly 
with  increasing  distance. 

The  contribution  of  this  cap  to  the  Vening-Meinesz  formula  is 


vc(5)  =  -  — - —  J*  Grad^St(5*t})  Ag(r?)  d  T  (rj) 

C 


(D.l) 


The  singularity  in  the  Vening-Meinesz  formula  stems  from  the  term  1/sin  in 
Stokes’  kernel  (B.24).  This  term  a’so  equals  2/£ ( §•  77) •  If  we  split  therefore 


st(s-”)=  »£5>  +  R<^> 


(D.2) 


then  the  contribution  of  R(§-rj)  leads  to  a  regular  integral 

vq( O  =  -  -±—  J*  Grad  R(5.77)Ag(Tj)dr<T7) 

1  4irG  “  * 


(D.3) 


The  contribution  of  2/j2(^.rj)  has  to  be  dealt  with  further: 


vi<!)  °  '  TTg"  <f  Giadf,  A8<”)dr(’’) 

C 


(D.4) 


If  this  is  done  in  a  satisfactory  way,  then 

VC(?>  =  VjU)  +  v2(5)  (D.5) 

W<_  need  now  some  vector  analytical  preparations.  Call  o(r,rj)  a  unit  vector 
which  is  tangential  to  f  in  the  point  %  and  is  contained  in  the  plane  spanned  by  | 
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and  77.  In  other  words,  a  (§,77)  is  orthogonal  to  5  and  co planar  with  §  and  77. 
Call  rd.rj)  the  vector  §x  a(5»7?).  We  view  a  (5, 77)  as  a  vector  which  is  im¬ 
bedded  in  3 -space.  Then  the  operators  Graci^  and  Lap^  can  be  applied  to  each  cf 
the  three  components .  The  following  can  be  verified  (cos  $  =  §.77 ) 

Grad  0(5,77)  =  —--7  t  (5 , 77)  t(5,  77 )T  (D.6) 

V  sin  ip 

LaP„  a(?,77)  =  -  - 1 —  cr(  5 , 77 )  (D.7) 

sinz  0 

We  do  not  give  a  detailed  derivation  of  these  formulas .  Chapter  6  in  Meissl 
( 1 971a) deals  more  extensively  with  the  underlying  concepts. 

Performing  the  differentiation  in  (D.  4)  gives 

1  n  COs3  4 

vi<5>  =  -  r~r  I  — y — “  «(5.i?)^(t?)dr(T,) 

2  v  G  c  s“*  $ 

Hence  by  (P.7) 

Vj(5)  =  ~~r  I  LaP-  ^(5,7?)  Ag(77)d  r  (77 )  + 

ZttG  ” 

C 

7  ,  1  -  COS^-y 

+  — -  f  - — — — 0(5,77)  Ag(7?)  d  r  (77) 

2nG  c  sin2  ip 

=  vxl(?)  +  vi2<?)*  say-  (°-8) 

Vj2(5)  is  already  regular.  The  integrand  is  even  bounded.  We  apply  Green’s 
second  formula  (component-wise)  toward  Vy(5)  and  obtain: 
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vu(?)  =  I  0(5,77)  Lap  ^(77)  dT  (77)  • 

^  IT  O 

c 

"  7T  #  a(?,T?)  (Grad  Ag(tj),  i/(T?))d3C(rj)  + 

^tr  o 

3C 

+  27q  #  Ag(T?)  (Grad^  a(  ?,  77),  y(r?) )  d  d  C  (rj) 
3  C 


The  last  integral  vanishes  because  of  (D.6)  and  the  orthogonality  of  r(|,  77)  and  1/(77). 
Summarizing  we  have  the  following: 

v  (5)  =  — L_  J  a(§,  77)  Lap  Ag(rf)  d  r  (77)  - 

2it  G 

C 


“  J  cr(5»  7?)  (Grad  Ag(77>,  y(i7>  )  d  3  0(77)  + 


2  ttG 


3  C 


i  (.  1  *  cos3  -|- 

—  J  - - *-  ct(| ,  77 )  Ag(77>  d  T  (77 )  - 


2ttG  1  Sin2  0 
C 

f  Grad^  R(  1*77)  Ag(7?)  d  r  (77) 
C 


(D.9) 


In  the  literature,  the  Vening-Meinesz  formula  is  usually  written  in  a  local 
ih  ,  a  coordinate  system.  In  this  notation  (D.9)  reads  as: 

2tt 

f 

2nG 


vc 


lp0  zn 

—  J  J  Lap  Ag(j^),  a)  sin  0  dtf  da 

”G  *=o  0=0  V sin  «  / 


2  IT 

-L_  sin  ij)n  f  (cos  <*\  3_Ag(0  it  a) 

170  0=0  \sin  a  ) 


4>  =  $l 


da  + 


^0  2W  o  .. 

1  (.  |»  1  -  cos^  $ 

rrG  *>  - ~ - 

t/)~ 0  a=0  sin  ip 

.  Vo  2 «■ 

J_  f  f  d  R(cos  ib) 

IT  G  *•  " 

lj)  =  0  0=0  dip 


(  cos  a\ 

Ag(0,  a)d$  do  - 

\sin  a  J 

(  cos  a\ 

Ag(tp,a)  sin  0  dji  do 

\  sin  o  / 


....  (D. 9a) 

In  the  main  portion  of  this  report  the  Vening-Meinesz  formula  occurs  in  (1.2X 
Ag(n)  has  then  to  be  replaced  by  Ag(tj)  plus  some  correction  terms. 
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Deflection  -  Formula 


»Ve  shall  he  concerned  with  t*e  following  term  contributing 
toward  the  deflection  formula  (1.2)  : 


4rR  G 


2  Grade  -*-1 -  Ag(q)  df  (q) 


rq.ri) 


Using 


Grad 


5  Y 
3  cos^-f 


fj  PIhT  sin4 V 


W— V  =  {— 

« i  sm  v  J  su 


—  O'  (],»)) 


...  (E.l) 


•  (E.2) 


fy  sin  y 


<y(j,^)  (E.3) 


we  may  write  the  cap  contribution  of  (E.l)  as  : 

qp(n  =  1  -g-  \  L(h(n)-h(?))2  Ag(a)]  Lap V  cr(|,7)j  dT(^) 

^  0  4rR  G  J  ‘(.sin  w  J 


[(»(-!)-»(!) )2Ag(7)]j 

- "V*  O’  (M)  d  r(fl) 

sin  y  J 

=  q1(j)  +  q2(f)  »  say* 


sin  y 


(E.4) 


q2(|)  is  already  regular.  Toward  qx(j)  Greens  second 
formula  will  be  applied  yielding 
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llff)  =  1  ““5^;  ^[(^(f))2^,)].^) 


“  ~  y  ~r~2^  <7(|^)  (  Grad,j[(h(^)-h(j))2  Ag(^), 

»  V  (^))  d  3C (i^) 


2>C 


-i  r  o  cos  y 

"TTTg  J  tOi^J-bCj))' ^Agty)]  — j-  <y(j,7)  d3c(7) 

3C  C  ^ 


...  (E.5) 


For  the  last  term  (D.6)  has  been  used.  The  applicability 
of  Greens  second  formula  has  to  be  justified  in  the  usual 
way  (letting  smother  small  cap  contract  toward  ^  ). 

The  first  term  in  (E.5)  contains  still  a  singular  integral. 

»7e  have: 

LaP^[(h(i2)-h(j)  )2  Ag(y)]  =  2  |  Grad  h(^)j2Ag(i£)  plus 

terms  containing  (h(^)-h(J))  as  a  factor.  (E.6) 

The  terms  containing  (h(r2)-h(J))  as  a  factor  cause  no 
singularity  problem.  However  the  term  2|Grad  h(^)f2  does. 

There  is,  however,  no  need  for  further  regularization. The 
term 

- i-n-  f  — ^-<5'(f»7)  J  Grad  h(n)|  2  Ag(n)  d  T(n)  (E.7) 

2rrG  £  sin^y  J  1  1  L 

can  be  combined  with 

-  J  Hrad  ^  .-y  j  Grad  h(q)  |  2  Ag(lf)  d  T(^)  (E.8) 

in  a  way  that  the  singularity  cancels.  (E.8)  is  contributing 
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to  the  first  term  of  the  deflection  formula  (1.2).  To  see 
this  consider  the  cap  contribution  of  (1.2),  split  St(j) 
according  to  (D.2)  and  note  the  definition  of  Gg(j)  in  (1.4). 
The  result  of  combining  (E.7)  and  (E.8)  is 


1 

— 

2  TC  it  Gr 


1 


1  -  co.5  V 
sin2  y 


<y(j»>?)  |  Grad  h(^)|^^g(^)  dP(i^) 


...  (E.9) 


which  is  now  regular. 

The  result  of  this  somewhat  involved  regularization  procedure 
is  summarized  in  section  3.5.4. 
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